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I .  INTRODUCTION 


In  many  simulation  applications,  it  is  required  to 
generate  bivariate  random  variables  which  have  identical 
marginal  distribution. 

For  example,  in  reliability  problems,  the  assumption  that 
two  components  have  independent  exponential  failure  times  is 
very  often  unrealistic,  and  much  effort  has  gone  into  deriving 
bivariate  exponential  random  variables  to  handle  these  situa¬ 
tions  [Gaver  (1972),  Olkin  &  Marshall  (1967),  Downton  (1970)]. 
Now  except  in  very  specific  physical  situations  it  may  be 
difficult  to  specify  the  complete  bivariate  distribution  of 
life  time  of  each  component.  However,  it  may  be  realistic 
to  specify  the  marginal  distributions  and  some  measure  of 
dependence  (usually  the  correlation  coefficient)  between  the 
life  time  of  each  component.  In  this  kind  of  situation,  we 
can  use  bivariate  random  vectors  having  given  marginal  dis¬ 
tribution  and  dependence  to  solve  the  problem  in  simulation. 

To  generate  these  vectors,  there  exists  some  previous  works 
and  existing  methods,  discussed  in  Section  II.  As  seen  in 
Section  II,  most  previous  work  is  specific  to  specified 
marginal  distributions  and  uses  inverse  transformation  methods 
as  a  basic  concept. 

An  example  is  the  recent  work  by  Johsnon  and  Tenenbein 
(1979),  to  generate  a  bivariate  random  vector  (X,Y)  which  has 
marginal  distribution  F^(x),  and  correlation  p ,  by  a 

weighted  linear  combination  method. 
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Define 


X  =  F^1(H1(U)) 

Y  =  F"1(H2(V)) 


or 


Y  =  F”1 (1  -  H2 (V) ) 

where  and  H2  are  the  cumulative  distribution  functions 
(c.d.f.)  of  U  and  V  respectively  and 


U  =  U' 

V  =  cU'  +  (1  -  c) V' 

where  U',  V'  are  i.i.d.  random  variables  with  probability 
density  functino  g(  ).  In  this  procedure,  F^,  F ^ ,  and  c 
are  specific  to  the  marginal  distribution  and  correlation 
desired.  The  functions  F^  and  F^  are  difficult  to  compute 
in  most  cases  and  the  weighting  factor  c  is  also  difficult  to 
calculate.  Moreover  most  of  the  work  in  univariate  random 
number  generation  has  been  aimed  at  avoiding  having  to  calcu¬ 
late  inverse  cumulative  distribution  functions  such  as 
F^(.)  and  F^(-).  These  are  the  reasons  why  many  proposed 
methods  are  specific  to  a  specified  marginal  distribution. 
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Again  special  properties  of  certain  random  variables  such 
as  infinite  divisibility  have  been  exploited  to  give  easily 
generated  bivariate  random  variables,  often  though  with 
limited  ranges  of  dependency.  One  very  clever  scheme  by 
Gaver  (1972)  to  generate  bivariate  exponential  random  varia¬ 
bles  uses  the  fact  that  the  sum  of  a  geometrically  distributed 
number  of  exponential  random  variables  (Y)  is  exponentially 
distributed  and  that  the  minimum  of  this  geometrically  dis¬ 
tributed  number  of  independent  logistic  random  variables  (Z) 
is  exponential.  Clearly  when  Y  is  large,  Z  is  small.  This 
scheme  is  of  course  very  specific  to  exponential  marginal 
distributions  and,  via  an  exponential  transformation,  to 
uniform  random  variables.  To  avoid  these  kinds  of  limitations 
and  to  make  the  generation  of  bivariate  random  variables 
simpler  and  more  automatic  in  simulations  we  develop  here  a 
scheme  presented  in  Jacobs  and  Lewis  (1977) . 

This  scheme,  the  mixture-truncation  method,  is  a  very 
general  tool  which  requires  only  that  a  method  be  available 
for  generating  random  variables  with  the  desired  marginal 
distribution.  The  mixture- truncation  method  scheme  for 
generating  bivariate  random  variables  is  as  follows. 

Let  F(x)  be  the  common  marginal  distribution,  of  the 
bivariate  random  variable  (Y,Z).  Let  p  be  the  desired 
correlation  between  Y  and  Z  (which  may  or  may  not  be  attain¬ 
able  generally  or  in  particular  with  the  mixture-truncation 
scheme) . 
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Define  the  transition  matrix  P  as 


P  = 


1  -  a. 


1  ~  a2  “2 


with  stationary  vector 


1  -  a. 


1  -  a. 


tr  P  =  tt  =  ( 


l-a^  +  l~a2,  l-a^  +  l-  o^' 


and  let  the  range  of  the  random  variable  X  with  distribution 
function  F(x)  be  x*  Then  generate  (Y,Z)  as  follows. 

1.  ( Initialisation) 

i)  Choose  an  "allowable"  xQ  from  (x.„,x  )  (the 

allowable  range  of  xQ  will  usually  be  smaller 
than  x  and  depends  on  p  and  F(x)). 
ii)  Set  =  F(xq)  ,  t\ 2  =  1  -  and  compute  and 

a  2  • 

iii)  Denote  by  X^  the  random  variable  X  truncated  to  the 

i 

left  of  xQ  and  by  X^  truncated  to  the  right  of 

x  . 

o 

2.  (Generate  Y)  .  Choose  Y  from  X^  with  probability  tt^ 
or  choose  Y  from  X2  with  probability 

3.  (Generate  Z) 

i)  If  Y  is  chosen  from  ,  choose  Z  from  X^  with 

probability  or  choose  Z  from  X2  with  probability 
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ii)  If  Y  is  chosen  from  X2,  choose  Z  from  X]L  with 
probability  1  -  or  choose  Z  from  X2  with 
probability  a2« 

In  Section  III,  we  will  show  that  the  correlation  p  between 
Y  and  Z,  if  xQ  is  properly  chosen,  is 


p 

E[[Y-E{Y}]  [Z  -  e<s>] 

°Y  °Z 

= 

6  M, 

where 


6  = 

-  (1  -  a2)  , 

M  = 

(  yl  “  2  77 1  17 2 

2 

°X 

P1 

=  E{X1),  u2  =  E{X2>  , 

2 

°X 

=  VAR  [X  ]  , 

Moreover,  the  generated  bivariate  random  vector  (Y, Z)  will 
have  marginal  distributions  F(x)  and  correlation  p  and  the 
joint  distribution  function  of  (Y,Z)  will  be 
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if 


F(y,z)  = 


F(z)  F  (y) 


z  <  x  ,  y  <  x 
—  o  —  o 


(1-0.) 

{o1  +  — - - [F  (z) -F  (xQ)  ]  }F  (y)  if  z  >  xQ,  y  <xq 


(1-a  ) 

(a.  +  - [F  (y)  -F  (x  )  ]  }F(z)  if  z<x  ,  y>x 

1  IT  2  °  “  O  c 


<*1*1  +  (1“°‘1)  [F  (z)  -F  (xq)]  —  if  Z  <  xq,  y  >  XQ 


+  { ( l-a„ )  + — (F(z)-F(x)  ] }  [F(y)-F(x  )  ] 
A  it  2  o  o 


These  relationships  will  be  developed  in  Section  III.  Note 
that  ( Y, Z )  may  be  continuous  or  discrete  random  variables 
or  mixtures  of  both,  though  in  this  thesis  we  concentrate  on 
continuous  cases. 

The  key  problem  in  this  very  simple  algorithm  comes  at 
the  initialization  steps  i)  and  ii) .  There  are  two  degrees 
of  freedom  in  the  selection  matrix  P  but  setting  =  F(xq) 
constraints  reduce  to  one  degree  of  freedom.  Specifying 
a  desired  correlation  further  constrains  the  degree  of  free¬ 
dom,  though  not  completely.  Subject  to  the  constraint  that 


and 

c*2  are  probabilities  there  may 

be 

a. 

no  values  xq  which  will  give 

(Y, Z)  with  correlation 

b. 

P 

one  value  x 

o 

c . 

a  range  x0- 
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The  main  numerical  problem  of  initialization  step  (i)  then 
is  to  compute  the  range  of  allowable  xq  for  a  given  p  and 
F  (x)  . 

The  main  statistical  problem  is  then  to  choose  which  of 

the  bivariate  random  vectors  (Y,Z)  indexed  by  xq  e  x0t0  use. 

An  alternate  solution  to  the  statistical  problem,  giving 

another  algorithm  is  then  to  let  xq  have  some  distribution  in 

the  range  Xq'  possibly  uniform  or  triangular,  this  not  only 

alleviates  the  problem  of  picking  a  particular  xq  e  xQ  ^ut 

it  also  smoothes  out  the  distribution  of  (Y,Z)  and  possibly 

makes  it  continuous.  The  computation  of  x  for  given  p  is 

o 

illustrated  in  subsequent  sections  for  uniform,  exponential 
and  gamma  marginals  for  (Y,Z). 

Before  doing  this  and  developing,  in  Section  III,  the 
results  already  given  here,  we  review  a  few  existing 
methods  for  generating  bivariate  random  variables  in  Section 
II.  These  are  methods  which  seem  fairly  tractable,  for  use 
on  a  computer. 
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II.  REVIEW  OF  EXISTING  METHODS 


There  are  several  existing  methods  for  generating 
bivariate  random  vectors,  but  most  of  them  are  specific  to 
particular  marginal  distribution  and  use  inverse  transfor¬ 
mation  method.  The  inverse  transformation  method  is  a  very 
useful  univariate  procedure  which,  unfortunately,  is  not 
possible  to  use  with  many  distributions  because  it  is  diffi¬ 
cult  and/or  uneconomic  to  compute  the  inverse  functions.  We 
survey  here  some  of  the  methods  which  are  germane  to  this 
thesis,  in  particular  concentrating  on  bivariate  random 
variables  with  uniform,  exponential  and  gamma  marginals.  First 
we  review  the  problem  of  determining  the  range  of  correlation 
coefficient  p  which  can  be  obtained  for  bivariate  distribution 
with  given  marginal  distributions,  again  considering  only 
the  continuous  case. 

A.  RANGE  OF  CORRELATION  COEFFICIENT  p 

Suppose  that  Y,  Z  are  random  variables  with  an  arbitrary 
joint  distribution  F(y,z)  with  finite  second  moments.  Then 
in  general  the  correlation  coefficient  p  can  take  all  values 
in  the  closed  interval  [-1,1].  But  with  specific  marginal 
distribution  F^(y)  and  F2(z),  the  class  of  all  F(y,z)  need 
not  attain  the  values  of  -1,  1,  of  p.  The  necessary  and 
sufficient  conditions  that  there  exist  determination  of  F(y,z) 
with  p  equal  to  1  and  -1  are  given  by  Moran  (1967)  as  follows. 
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i)  there  exist  constants  a  and  6  such  that  aY  +  3  has 
the  same  distribution  as  Z, 
ii)  the  distribution  of  Z  is  symmetrical  about  its 
mean. 

To  see  this,  rescale  Y  to  have  the  same  mean  and  variance 

as  Z.  If  there  exists  an  F(y,z)  such  that  p  =  1  we  have 
2 

E[Z-Y]  =0,  so  that  Y  =  Z  with  probability  one.  Then  Y 

must  have  the  same  distribution  as  Z.  On  the  other  hand 

if  there  exists  an  F(y,z)  such  that  p  =  -1  we  shall  have 
2 

E[Z  +  Y]  =  0,  and  Y  =  -Z  with  probability  one.  Thus  if  both 
bounds  are  attainable  Z  must  have  a  symmetric  distribution. 

Given  general  marginal  distributions  F^(y)  and  F2(z), 
Mardia  (1970)  showed  Frechet  bounds  as 

max[0,F1  (y) +F2  (z) -1]  <_F(y,z)  <_  min  [F1  (y)  ,  F2  (z)  ]  (II-A-1) 

From  this  we  can  find  the  range  of  possible  values  of  p. 

For  simplicity  we  now  confine  our  consideration  to  distribu¬ 
tions  F(y,z)  of  positive  random  variables  whose  deriva¬ 
tives  F|(y)  and  F2(z)  are  strictly  positive  for  y  >  0,  z  >  0, 
respectively.  Suppose  also  that  the  variates  are  scaled  to 
have  unit  variances.  Let  (u) ,  G2 (v)  be  the  inverse  func¬ 
tions  of  F^ (y) ,  F2(z),  i.e. 

F i  [Gj_  (w)  ]  =  w 

where  0  <  w  <  1,  i  =  1,  2.  Then  the  correlation  coefficient 
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between  Y  and  Z  is  given  by  the  equation 


00 


00 


p 


J  /  y  Z  dF(y,z)  -  E[Y]  E[Z] 
0  0 


( II-A- 2a) 


1  1 

/  /  G,  (u)  G  (v)  dK(u,v)  -  E[Y]  E[Z] 

0  0  X 


(II-A- 2b) 


where  K(u,v)  is  the  joint  distribution  of  the  quantities 
U  =  F^(y),  V  =  F 2(z) •  Then  U,V  are  jointly  distributed  on 
the  square  0  _<  U  £  1,  0  <_  V  <_  1,  in  such  a  way  that  the 
marginal  distributions  are  uniform  on  the  unit  intervals. 
From  expression  (II-A-1)  the  minimum  correlation  is  attained 
when  the  probability  is  concentrated  uniformly  on  the  line 
U  +  V  =  1.  The  minimum  value  of  p  then  is 


1 

J  G1(u)  G2(1-u)  du  -  E[Z]  E  (Y  ]  (II-A-3) 


pmin 


0 


The  corresponding  F(y,z)  will  be  a  singular  distribution  with 
all  the  probability  concentrated  on  the  line  F^(y)  +  F  ^{t.)  =  1. 
In  fact  Y  and  Z  are  an  antithetic  pair.  The  maximum  value 
of  p  is  attained  when  the  probability  is  concentrated  uni¬ 
formly  on  the  line  U  =  V.  Then  the  maximum  value  of  p  is 


P 


max 


0 


1 

J  Gx(u)  G2(u)  du  -  E  ( Z  ]  E[Y] 


(II-A-4) 
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with  the  corresponding  probability  concentrated  on  the  line 
F1 (y)  =  F2  (z) * 

By  using  this  result  we  will  get  lower  and  upper  bounds 
of  the  correlation  coefficient  for  uniform,  exponential  and 
gamma  marginal  cases.  In  the  uniform  marginal  distribution 
case,  we  can  see  p  .  =  -1  and  p  =1  from  Moran's  condition 

ii) ,  i.e.,  uniform  distribution  is  symmetric  about  its  mean. 
For  the  exponential  marginal  distribution  case  suppose  Y  and 
Z  have  exponential  distributions  with  unit  mean  and  variance. 
Then 


F1(y)  =  1  -  e"y,  F2(z)  =  1  -  e  2 

and  the  inverse  functions  are 

( u)  =  -  In  (1  -  u)  , 

G2(v)  =  -  In  ( 1  -  v)  . 

From  equations  (II-A-3)  and  (II-A-4) 

1 

pmin  =  J  Gl(x)  G2(1_x)  dx  "  E[Y]  E[Z] 

1  2 
=  f  In  x  In  ( 1-x)  dx  -  1  =  1  - 

0  6 

;  -0.64493 
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1 

p  =  f  G,  (x)  G„(x)  dx  -  E[Y]  E[Z] 

max  „ ‘  1  i 


1 

/  In  x  In  x  dx  -  1  =  1 

0 


The  p  .  can  be  attained  when  all  Y  and  Z  are  concentrated 
min 

on  the  line  e  2  +  e  =1.  Also  the  p  can  be  attained 

max 

when  all  Y  and  Z  are  concentrated  on  the  line  e  2  =  e  , 


For  the  gamma  marginal  distribution  case,  suppose  that 
Y  and  Z  have  gamma  type  distributions  with  probability  density 
functions 


fi(y)  ■  rToT  e"Y  2“'1  “  "  0  •  y  i  0 

f2(z)  -  rwr  e’z  28-1  6  =■  0  -  z  1 0  • 


Then 


E[Y]  =  a,  E  [Z  ]  =  6,  VAR  [Y  ]  =  a  VAR[Z]  =  6  . 


The  Pmj_n  will  be  attained  when  the  probability  density  is 
concentrated  on  the  line 


1 

r  (a) 


0 


o-l 

x 


dx  + 


r  (e) 


-x 


.6-1 


dx 


1  . 
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This  defines  y  uniquely  as  a  function  of  z  which  can  be 
written  y  =  A(z)  .  The  Pmj_n  is  then,  on  rescaling  the 
covariance 


p 

min 


{  /  u  A (u)  f,  (u)  du  -  ctB)/(aB)  1//2 
0 


When  a,  8  become  large,  p  .  tends  to  -1.  And  the  p 

min  max 

will  be  attained  when  the  probability  density  is  concentrated 
on  the  line 


1 

r  (a) 


( y  -X  a-1  , 

J  e  x  dx 

0^ 


1 

r  (B) 


2  i 

r  -x  a-1  , 

J  e  x  dx 

0 


This  defines  y  as  a  function  of  z  which  can  be  written 
y  =  B(z).  The  p  then  is,  on  rescaling  the  covariance 

IT13.X 


p 

max 


{  {  u  B  (u)  f ,  (u)  du  -  a8}/(aB)1//2 
0 


when  a  =  B,  p  tends  to  1.  Schmeiser  and  Ram  Lai  (1979) 
max 

showed  the  obtainable  correlations  between  random  variables 
Y  and  Z  having  gamma  marginal  distribution  with  density  function 

a  •  -1 

fj_(x)  =  [  (x/  8^)  1  exp  (-x/B  j_)  /  (8  j_r  (a^)  ] 

for  X  >  0,  >  0,  Bi  >  o,  i  =  1,  2. 
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The  Figures  (Il-a)  and  (Il-b)  show  the  obtainable 
correlation  as  a  function  of  a given  =  1  and  5,  for 

Bi  ■ 


I 


a* 


Figure  Il-a.  Obtainable  correlation  as  a  function  of 
a ^  for  a1  =  1,  81  =  1,  by  Schmeiser 
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a,=i 


0  5  30  oo 

Oil 


Figure  Il-b.  Obtainable  correlation  as  a  function  of  a2 
for  =  5,  6^  =  1,  by  Schmeiser 
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B.  GENERAL  METHOD  REVIEW 


1 .  Johnson  and  Tenenbein's  General  Method 

A  general  method  of  constructing  a  bivariate  dis¬ 
tribution,  whose  marginal  distribution  functions  are  F^(x) 
and  F2(y),  is  proposed  by  Nataf  (1962),  can  be  represented 
as  follows. 

General  Method 

i)  Consider  any  two  continuous  random  variables  U 
and  V  with  probability  density  functin  h(u,v). 

ii)  Let  X'  =  H^(u)  and  Y'  =  H2 (v) ,  where  H^(u)  and 
H2(v)  are  the  cumulative  distribution  functions 
of  U  and  V,  respectively, 

iii)  Define 


X  =  F11(X’)  =  F"1(H1(u)] 


(II-B-1) 


and 


Y  =  F^1  ( Y  ’ )  =  F21[H2(v)] 


(II-B-2a) 


or 


Y  =  F^1  ( 1  -  Y '  )  =  F21(1-H2(v)]  (II-B-2b) 
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Since  X',  Y'  and  1-Y'  are  uniformly  distributed  over  the 
interval  [0,1],  X  defined  by  expression  (II-B-1)  and  Y 
defined  either  by  (II-B-2a)  or  (II-B-2b)  will  have  a  joint 
distribution  whose  marginal  distribution  functions  are  F^(x) 
and  F^Cy).  The  method  is  probably  what  one  would  think  of 
first  but  again  involves  inverses.  Based  on  this  general 
method,  Johnson  and  Tenenbein  (1979)  developed  two  proce¬ 
dures  for  generating  (and  also  constructing)  bivariate  dis¬ 
tributions  whose  marginal  distributions  and  measures  of 
dependence,  as  given  by  Kendall's  T  or  the  grade  correlation 
coefficient  pg,  are  specified.  Note  that  these  Kendall's  T 
and  grade  correlation  coefficient  pg  are  not  the  same  as  the 
ordinary  product  moment  correlation  coefficient  p  we  have 
been  considering;  these  measures  are  discussed  by  Kendall 
(1962)  and  Kruskal  (1958)  and  can  be  defined  as  follows.  Let 
X  and  Y  be  continuous  random  variables  having  some  joint 
probability  density  function.  Let  (X^,Y^),  ^2^2)  and 

(X^rY^)  be  three  independent  pairs  of  observations  having 
the  same  joint  density  function,  then 

T  =  2  P[(X1-X2) (Y1-Y2)  >  0]  -  1 

Ps  =  6  P [ (X1-X2) (Y1-Y3)  >  0]  -  3 

The  first  procedure  for  generating  bivariate  pairs,  called 
the  WLC  (weighted  linear  combination) ,  defines 
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u 


U' 


(II-B-3) 


and 


V 


cU'  +  (l-c)V 


( II-B-4 ) 


where  0  <_  c  <_  1,  U'  and  V'  are  independent  and  identically 
distributed  random  variables  with  probability  density  func¬ 
tion  g(»).  Then  X  and  Y  are  obtained  from  the  general 
method  using  equation  (II-B-1)  and  (II-B-2a)  if  the  depen¬ 
dence  measure  is  positive  or  equations  (II-B-1)  and  (II-B-2b) 
if  the  dependence  measure  is  negative. 

In  order  to  apply  the  WLC  procedure,  we  must  obtain 
expressions  for  H^(u),  H^Cv),  as(u,v)  and  T(u,v),  in  terms 
of  c  and  g(.)  .  The  values  of  H-^(u)  and  (v)  allow  us  to 
apply  the  general  method  for  a  given  choice  of  c  and  g(t). 

The  expressions  for  T (u,v)  and  as (u,v)  allow  us  to  specify 
c  for  a  given  choice  of  g(-)  in  terms  of  the  required  value 
of  either  T  or  c  .  From  equations  (II-B-3)  and  (II-B-4) 


u 


H]_(u) 


j  g  ( t)  dt 


(II-B-5) 


—  CO 


H2(v) 


/  /  g(u')  g(v')  du'dv'  (II-B-6) 
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where 


=  {  (u ' , v ' )  :  cu '  +  ( 1-c) v '  <_  v) 


and  the  joint  density  function  of  u  and  v  is 


h  (u,  v) 


1  ,  .  ,v  -  cu, 

j-rj  g(u) 


(II-B-7) 


OO  CO 


(u,v)  =12  j  J  (u) ^ (v) h (u, v) dudv  -  3  (II-B-8) 


—  00  —  00 


T(u,v)  =  4  }  G_  t)  g  (t)  dt  ,  (II-B-9) 


where 


g2(t)  =  /  g(w  +  t)  g(w)  dw 

—  00 

00 

G2(t)  =  /  g2(x)  dx 

—  co 

Using  equations  (II-B-5) ,  (II-B-6) ,  (II-B-7) ,  (II-B-8)  and 

(II-B-9),  we  have  to  evaluate  H1(u),  H2 (v) ,  h(u,v),  T(u,v) , 
and  pg(u.v)  for  all  values  of  c  and  for  specified  g(*)* 

As  Johnson  and  Tenenbein  mentioned  in  their  report,  these 
integrals  are  quite  tedious  to  perform,  so  they  gave  some 
computation  results  in  their  report. 
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The  second  procedure,  called  TVR  (trivariate  reduction) 
is  discussed  by  Mardia  (1970) .  In  this  case  U  and  V  are 
defined  as 

U  =  U'  +  BZ'  (II-B-10) 

V  =  V'  +  BZ'  (II-B-11) 

for  0  <_  3  <  00 ,  and  U',  V',  and  Z'  are  independent  and  iden¬ 
tically  distributed  random  variables  with  probability  density 
function  g ( • ) . 

Like  the  WLC  procedure,  one  needs  to  define  H^(u), 

H„ (v) ,  h(u,v),  p  (u,v)  and  T(u,v)  as  a  function  of  3  and  g(*). 

^  s 

From  equations  (II-B-10)  and  (II-B-11)  it  follows  that 

H1(u)  =  h2(v)  =  /  Jg(u')  g(z')  du'  dz '  (II-B-12) 


where 


=  {(u',z'):  u'  +  3z'  <_u}. 


00 

h(u,v)  =  /  g (u  -  3z) g (v  -  3z)  g(z)  dz  (II-B-13) 

—  00 


and 
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T  (u,  v) 


( II-B-14 ) 


00 

=  4  /  [G2(et) ]2  g2(t)  dt  -  1 

—  00 


where 


g2  ( t)  =  /  g(w+t)  g(w)  dw 

—  oo 

oo 

G2(t)  =  /  g2 (x)  dx 

—  00 

using  equations  (II-B-12),  (II-B-13),  (II-B-14)  and  (II-B-8) 

the  same  computations  are  needed. 

In  this  general  WLC  and  TVR  methods,  there  are  two 
problems.  One  is  the  need  for  inverse  transformations.  The 
other  is  the  need  to  generate  coupled  random  variables  to 
create  dependence  after  the  inverse  transformation  is  applied. 
Doing  this  by  linear  combinations  is  not  necessarily  the 
simplest  and  most  felicitous  with  respect  to  calculation  of 
correlations  and  range  of  correlations. 

A  significant  simplification  can  be  achieved  by 
obtaining  a  "smooth"  bivariate  uniform  pair  (U,V)  whose 
correlation  can  be  varied  between  the  limit  of  correlation 
which  can  be  obtained  for  uniform  marginals.  These  limits 
are  (-1,1)  since  the  antithetic  pair  (U,l-U)  has  correlation 
-1. 
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Note  that  the  problem  is  essentially  one  of  obtaining 
a  bivariate  uniform  random  variable  with  positive  correlation 
because 


corr (U, V)  =  p 

u,v 


~PU,1-V  =  corr ( U, 1-V) 

“P1-U,V  =  corr (1-U, V) 


Two  fairly  simple  schemes  for  obtaining  bivariate  uniform 
distributions  are  discussed  in  the  next  section. 

C.  REVIEW  OF  SOME  SPECIFIC  METHODS 

The  mixture-truncation  method  is  a  general  algorithm 
for  generating  bivariate  pairs  (Y,Z)  with  given  marginal 
distributions  and  requires  only  the  availability  of  a  generator 
for  the  marginal  distribution  and  calculation  of  constants. 
Unless  randomization  is  used,  however,  it  does  give  a  bi¬ 
variate  distribution  with  a  discontinuous  joint  distribution, 
as  in  Section  (III-D) .  This  could  be  a  disadvantage.  It 
is  interesting  therefore  to  compare  it  to  other  generations 
which  are  specific  to  the  given  marginal  distribution  and 
we  do  this  here  for  the  cases  for  which  the  mixture-truncation 
method  is  illustrated  in  this  thesis,  namely  uniform,  exponen¬ 
tial  and  gamma. 

There  are  many  schemes  available  for  uniform  and  exponen¬ 
tial  cases  though  it  should  be  noted  that  smooth,  simply 
generated  pairs  with  easily  calculated  correlations  have 
only  recently  been  available.  In  particular  we  give  here  two 
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very  recent  schemes  for  smooth  bivariate  uniforms  which  are 
easily  generated  and  whose  correlations  can  be  calculated; 
similarly  two  schemes  for  exponential  which  are  essentially 
the  uniform  schemes  after  transformation.  When  one  comes 
to  gamma  marginals  one  is  again  in  difficulty  and  we  describe 
a  recent  scheme  by  Schmeiser  and  Ram  Lai  which  involves 
however,  computation  of  the  inverse  gamma  distribution 
function  and  very  difficult  initialization  computations. 

The  mixture-truncation  method  is  quite  competitive  here.  In 
fact  if  one  can  compute  the  inverse  gamma  distribution  then 
bivariate  gamma's  can  be  computed  as  [Y  =  F  "'"(U),  Z  =  F  (V)  ] 
where  U  and  V  are  a  bivariate  uniform  generated  by  the  schemes 
mentioned  above;  this  would  be  simpler  for  any  marginal  dis¬ 
tributions  than  some  of  the  general  schemes  mentioned  in  the 
previous  sections.  Comparisons  of  the  generating  schemes  given 
in  this  section  with  the  mixture- truncation  method  are  given 
in  Section  IV-D  and  Section  V-D. 

1 .  Lawrance  and  Lewis's  Bivariate  Uniform 

Lawrance  and  Lewis  (1979)  show  that  a  simple  trans¬ 
formation  of  the  NEAR(l)  (New  Exponential  Autoregressive) 
process  gives  a  two-parameter  family  of  Markovian  random 
variables  with  uniform  marginal  distributions.  This  method 
generates  a  correlated  uniform  pair  as  a  multiplicative  mix¬ 
ture  of  uniform  variables,  where  the  parameters  a  and  6  take 
on  values  between  zero  and  one,  with  the  condition  that  they 
not  both  be  one.  Let  Y  be  a  uniform  [0,1]  random  variable, 
and  define 
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wp  a 


£  Y 


3 


Z 


wp  1-a 


where 


U 


wp  1-3 

1  -  (1-a) 3 


U 


(1-a) 3 


where  U  is  an  independent  identically  distributed  uniform 
(0,1)  random  variable.  The  correlation  structure  is  defined 
as  follows  as  a  function  of  a  and  3; 


3  ,  a  3 _ 

2  +  3  1  +  ( 1— a )  3 


P 


Y,Z 


The  model  can  be  reduced  to  a  one-parameter  model  by  suitably 
relating  a  and  3,  e.g.,  a  =  3  and  all  positive  values  of  p 
can  be  obtained.  Consequently  Y  and  1-Z  will  have  a  full 
range  of  negative  correlations.  A  generating  procedure  for 
this  bivariate  pair  of  uniform  random  variables  is  as  follows. 

Generating  Procedure 

1.  (Initialization) .  For  given  correlation,  find  suitable 
parameter  values  a  and  3 . 


2.  Generate  U.  ,  set  Y  •*-  U.  ;  r  =  1-a,  P  «*  ^ — -t-t 

-L  _L  _L  v  -L  ot ;  p 

3.  Generate  U2» 
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U2 

4.  If  U2  <  r,  set  U  «-  —  and  go  to  8. 


5.  Otherwise  U  *■ 


U2-r 
1  -  r* 


6.  If  U  <_  P,  set  Z  -jTp  an<^  9°  to  10* 

7.  Otherwise,  set  Z  Y6  ( ra  and  go  to  10. 


8.  If  0  <  P,  set  Z  +■  ^  and  go  to  10. 

9.  Otherwise,  set  Z  +■  ( • 

10.  Go  to  2  and  repeat  until  the  desired  number  of  bi¬ 
variate  uniform  variables  (Y,Z)  are  obtained. 

The  program  is  listed  in  the  appendix  and  the  scatter  plot 
and  bivariate  histogram  of  resulting  output  are  in  Section 


IV-D. 


2 .  Bivariate  Uniform  By  Transformation  of  Gaver's 

Bivariate  Exponential 

In  Gaver's  bivariate  exponential  to  be  discussed  in 
Section  II-C-4,  we  can  define  a  bivariate  exponential  random 
vector,  ( Y , Z )  with  correlation  p  =  -  P/2  as  follows: 


Y 


In  [ 


u1//N  (1-P) 


1-P 


(II-C-1) 


and  Z  is,  conditional  upon  N  =  n,  a  gamma  random  variable 
with  shape  parameter  n  and  mean  n(l-p),  where  N  is  a  geometric 
random  variables  with  probability  density  function 

X”  1 

f(x)  =  p(l-p)  ,  x  =  1,  2,  ...  and  U  is  a  uniform  random 

variable  over  [0,1].  From  this  we  can  generate  a  pair  of 
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uniform  random  variables  (V,W)  by  transformation.  Since  we 
know  that  to  generate  an  exponential  random  variable  X  from 
a  uniform  random  variable  U[0,1),  we  use  the  following 
transformation: 


X  =  -  In  ( 1  -  U  ) 


(II-C-2) 


The  resulting  X  has  the  exponential  distribution  with  unit 
mean.  From  equations  (II-C-1)  and  (II-C-2) ,  we  can  generate 
V,  a  uniform  random  variable  over  [0,1]: 


V 


=  (1-p)  U 

1  -  P  U 


1/N 

1/N 


Further  also  we  can  generate  W,  a  uniform  random  variable 
over  [0,1]  : 


1-P 

W  =  (  n  U. ) 

.  .  l 

1=1 

This  follows  since  we  know  that  Z  has,  conditionally  upon 
N  =  n,  the  gamma  distribution  with  shape  parameter  n  and 
mean  is  n(l-p) . 

In  general  the  resulting  V  and  W  will  be  negatively 
correlated  and  (V,l-W)  will  be  a  positively  correlated  pair. 
Because  the  correlation  structure  will  be  changed  by  trans¬ 
formation,  the  resulting  (V,W)  need  not  have  the  same 
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correlation  with  (Y,Z).  The  correlation  coefficient  of 
(V,W)  will  be  computed  as  a  function  of  P  as  follows. 


p  =  12  E[VW]  -  3 


where 


E[VW]  =  en  [E  [VW  |N  =  n]  ] 


ilrPJLP. 

1  -  P  u 


1/n 

l7n 


du 


[  /1u1-P  du ] n ] 
0 


Unfortunately  this  computation  of  p  as  a  function  of  P  is 
difficult.  As  an  example  we  used  here  the  same  function  for 
p  as  holds  in  the  exponential  case. 


Generating  Procedure 
1.  (Initialization) 

i)  Compute  P  for  given  correlation  p 
ii)  Choose  N  from  a  geometric  distribution  with 
parameter  1-P 


2. 


Generate  a  uniform  [0,1]  random  variable  and  define 


(1-P)  U. 


1/n 


V  = 


1  -  P  U. 


1/n 


3.  Generate  N  =  n  uniform  [0,1]  random  variables  U^, 
i  =  1,  ...,  n,  and  define 
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w  =  (  n  u. ) 1-p 

i-1  1 

4.  Deliver  (V,W)  and  go  to  2  until  a  sufficient  number 
of  bivariate  pairs  have  been  generated. 

The  program  is  listed  in  the  appendix  and  the  scatter  plot 
and  bivariate  histogram  of  resulting  random  vectors  are  in 
Section  IV-D. 

3 .  Marshall  and  Olkin's  Positively  Correlated 

Bivariate  Exponential 

Suppose  X  is  the  age,  or  length  of  service  of  the 
first  device  at  the  time  of  death,  governed  by  two  indepen¬ 
dent  Poisson  processes  Z-^(t),  Z^2(t)  with  parameters  A^, 

A ^2  respectively  and  Y  the  same  of  the  second  device, 
governed  by  two  independent  Poisson  processes  Z2(t),  Z^2(t) 
with  parameters  A 2,  A^2,  resPectively •  Further  suppose  the 
second  device  is  placed  in  operating  after  a  time  5  later 
than  the  first  device.  Then  the  joint  distribution  of  (X,Y) 
is  defined  as 

P  [X  >  x,  Y  >  y]  e  F  (x,  y) 

■  exp{-A-^x-A2y  -  A-^max  [x,y+min  (x,  6)  ]  }  if  6  _>  0 

-  exp{-A^x- A2y-A12max [y ,x+min (y , 6 ) ] }  if  6  <_  0 

( II-C-3) 

If  <5  =  0,  then  equation  (II-C-3)  reduces  to 
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P  [X  >  x, Y  >  y] 


F  (x ,  y ) 


=  exp{-A^x  -  X^y  ~  *12  max(x/YH 

Note  that  since  the  marginal  distributions  do  not  depend  on 

6, 


E[X]  = 


Xl+X12  ' 


VAR [X]  = 


(*1+  *12) 


2  ' 


E  [Y]  = 


*2  +*i2' 


VAR  [Y  ]  = 


( A2  +  X12^ 


We  also  have,  for  6  >  0 


E[XY]  = 


X12  e 


^A1+X12^A2  +  X12^  ^X1  +  X12^X2  +  X12^ 


and  the  correlation 


corr(X,Y)  =  p 


X/Y 


X12  ~  ^Xl+X12^  ^ 


(II-C-4a) 


where 


X  X1  +  X2  +  X12 ' 
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If  6  =  0,  then  the  equation  (II-C-4a)  reduces  to 


P 


xy 


(II-C-4b) 


This  characterization  can  be  represented  in  terms  of 
independent  random  variables.  For  6  >_  0 ,  if  there  exist 
independent  exponential  random  variables  U^,  U^r  and 
with  respective  parameters  X X2 1  ^^.2'  ^12'  t^ien 


X  =  mindJ^U^) 


min(U2,U3-6) 

if 

U3  i  6 

-  min(U2,U4) 

if 

U3  <  6 

It  can  be  verified  formally  from  the  relation 


P[X>x,Y>y]  =  P  [U^  >  x,  U2  >  y  ]  (P  [U^  >  x,U3  >  y+6  |  >  6  ] 

x  P[U3  >  6]}  +  P[U4  >  y]P[U3  >  x|U3  <  6]P[U3  >  6] } 

In  the  6=0  case,  the  representation  yields 


X  =  min(U1,U3) 
Y  =  min(U2,U3) 
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This  bivariate  distribution  has  a  line-discontinuity,  if 
E  [X]  =  E  [Y]  ,  along  the  line  x  =  y,  i.e.,  in  the  case  where 
x  =  y  =  •  Thus  this  bivariate  exponential  is  not  as  smooth 

as  others.  In  particular  the  NEAR(l)  process  of  Lawrance 
and  Lewis  generates  bivariate  exponentials  with  a  continuous 
density  function.  But  if  we  consider  the  case  6  >_  0 ,  the 
line-discontinuity  will  be  removed.  In  either  case,  the 
resulting  (X,Y)  pairs  always  have  positive  correlation  as 
shown  in  equations  (II-C-4a)  and  (II-C-4b) .  If  we  use  these 
methods  to  generate  positively  correlated  random  vectors 
(X,Y)  for  exponential  marginal  distribution  with  unit  mean 
and  given  p,  we  have  to  compute  A^ ,  A ^  and  A^  f°r  9iven 
correlations  to  generate  independent  exponential  random 
variables.  For  simplicity  we  will  consider  the  method  with 
6  =  0.  From  the  relationship  A  =  A-^  +  A£  +  A^2' 


E[X] 


*1  + 


l12 


1 


E[Y] 


a2  + 


12 


P 


1 


we  can  compute  A^,  an<^  as  functi°ns  °f  P* 

A 12  =  2p/  (1 +  p ) 


A 


1 


1 


A 


12  ’ 
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Generating  Procedure 


1.  (Initialization) .  Compute  A^ ,  A^  and  A^2  for  given 

0  <_  p  <_  1 . 

2.  Generate  three  independent  exponential  random  variables 

,  E2  and  with  unit  mean. 

3.  Rescale  the  independent  random  variables  and  set  as 


U1  "  El/Al 


u2  -  e2/a2 


°3  "■  E3/X12 


4.  Define  X  and  Y  as 


X  =  min[UL,U3] 

Y  =  min[U2,U3] 

5.  Deliver  (X,Y)  and  go  to  2  until  a  sufficient  number 
of  bivariate  pairs  have  been  generated. 

The  program  is  listed  in  the  appendix  and  the  scatter  plot 
and  bivariate  histogram  of  resulting  random  vectors  are 
shown  in  Section  V-D. 
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4 .  Gaver 1 s  Negatively  Correlated  Bivariate  Exponential 

This  negatively  correlated  bivariate  exponential 
is  generated  from  the  following  considerations.  Suppose  a 
particular  system  has  N  defective  elements;  here  we  assume 
N  has  geometric  distribution  with  probability  density 
function 


x-1 


f (x)  =  p(l  -  p) 


x  =  1.  2 


The  time  to  failure  of  each  defective  element  is  T.  with 


l 


density  function  F(t).  Further  suppose  repair  is  carried 
out  after  all  defective  elements  are  discovered.  If  we 


define  R  as  repairing  time  which  is  distributed  as  Gn(x) 


provided  N  =  n,  and  define  T  as  minimum  of  T^,  then  the 
joint  distribution  function  of  T  and  R  is  obtainable  from 


P  [T  >  t ,  R  <_  x]  =  l  (1  -  F(t)  ) 

n=l 


where 


(1  -  p)  P 


n-1 


0<p<l,  n  =  1,  2 


P 


n 


Then 


E[R|T  >  t] 


E  [R] 


1  -  P  [1  -  F(t)  ] 
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where  E[R]  is  the  expected  value  of  an  element  repair  time. 
Moreover, 


E  [R  |  T  =  t] 


E  [R]  [1  +  Q-~p-(  1  -  <|>(t)  )  ] 


where 


*(t) 


P [T  <  t] 


F(t) 


1  -  P[1  -  F  ( t )  ] 


From  this  regression,  we  know  that  if  T  is  long,  R  is  short 
and  vice  versa;  the  negative  relationship  is  clearly  present. 
Further  we  can  calculate  Cov[R,T)  as 


P  E[R]  E  [T]  (1  -  3 


Cov [R, T) 


0 


< 


where  T(l)  has  the  distribution  of  the  smallest  of  a  sample 
of  two  from  4>  ( t )  . 


It  will  not  be  shown  that  one  can  select  F(t)  =  l-F(t) 


in  such  a  way  as  to  induce  exponentially  distributed  time 
to  system  failure,  T.  Since  we  know  that  exponential  G, 
geometrically  compounded,  yields  exponential  R,  the  outcome 
will  be  a  (T,R)  pair  with  exponential  marginals  and  negative 
correlation. 


40 


If 


F  ( t ) 


1  -  F  ( t ) 


satisfies 


l  [F  ( t)  ]  n  pn_1  (1-p) 


e 


-xt 


n=l 


then  T  is  exponential  with  mean  1/A.  Solution  for  F 
yields , 


1 


F  (t) 


t  >  0 


p  +  (1-p)  eAt 


which  is  a  logistic  distribution,  left  truncated  at  t  =  0. 
If  X  =  1,  then  T  is  exponential  with  unit  mean,  while 
E  [T^  ]  =  1/2  ,  so 

corr  (R,  T)  =  -  2.  E  [R]  E  [T  ]  . 

If  G  is  chosen  so  as  to  make  R  exponential  then  it  may  be 
shown  that 


corr (R,T) 


Consequently,  a  greatest  lower  bound  for  the  correlation  in 


this  model  is  -  -j. 
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Generating  Procedure 

1.  If  0  <  p  <  -0.5,  Set  P  -e  - 2p 

2.  Generate  a  geometric  random  variable  N  with  parameter  1-P 

3.  Generate  a  uniform  (0,1)  random  variable  U  and  define 


Y  as 


Y 


4.  Generate  a  gamma  random  variable,  Z,  with  shape 
parameter  n,  and  mean  is  n(l-p) 

5.  Deliver  (Y,Z)  and  go  to  2  until  a  sufficient  number 
of  pairs  are  obtained. 

For  obtaining  negatively  correlated  bivariate  exponential 
random  vectors,  this  algorithm  is  very  simple  and  one  of  the 
few  available.  Moreover  the  correlation  is  known  explicitly. 
The  program  is  listed  in  the  appendix  and  scatter  plot  and 
bivariate  histogram  of  resulting  random  vectors  are  shown 
in  Section  V-D. 

5.  Arnold's  Vibariate  Gamma  Generator 


Arnold  (1967)  developed  a  Trivariate  reduction  method 


to  generate  bivariate  gamma  random  vectors  having  positive 
correlation  with 
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where  and  a.^  are  the  given  shape  parameters  for  the  marginal 
distributions.  Letting  "gamma  (a,g)"  denote  the  gamma  dis¬ 
tribution  with  shape  parameter  a  and  scale  parameter  g , 
so  that 


f(x)  =  FrT^T  ‘f’”'1  exp('x/B) ' 

where  r(a)  is  the  gamma  function  and 

E[x]  =  a  g  .  VAR [ x ]  =  ag2, 


the  trivariate  reduction  algorithm  proceeds  as  follows  for 
given  shape  parameter  values  a^,  a2  and  p  to  generate  unit 
scale  gamma  variate  i.e.,  g  =  1. 

Generating  Procedure 


1. 

Generate 

G1  ~ 

gamma 

(ct^  -  P  ( 

a/2  ,, 

a1a2)  /l) 

2. 

Generate 

G2  - 

gamma 

(a2  -  p  ( 

a/2 

a1a2)  ,1) 

3. 

Generate 

G3  ~ 

gamma 

(p  ( a]_a2 

)1/2,d 

4. 

Define  Y 

and 

Z 

Y 


G1  +  G3 


Z 


G2  +  G3 
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These  Y  and  Z  are  unit  gamma  variates  with  6=1,  but  Y  and 
Z  can  be  multiplied  by  6-^  and  &2'  respectively,  to  obtain 
any  desired  scaling.  In  the  algorithm,  is  a  common  com¬ 
ponent  of  both  Y  and  Z,  thus  inducing  positive  correlation. 

When  a  positive,  but  not  extreme,  correlation  is 
needed,  trivariate  reduction  is  an  excellent  algorithm, 
using  univariate  gamma  generators.  If  the  marginal  distri¬ 
butions  have  the  same  shape  parameter  value  a,  then  the 
correlation  limitation  is  0  <_  p  <  1.  But  in  the  other  case 
the  upper  limit  is  smaller  than  1. 

Note  that  this  algorithm  exploits  the  infinite 
divisibility  of  the  marginal  gamma  variables  and  is  applicable 
to  any  infinitely  divisible  marginal.  It  is  commonly  used, 
for  instance,  to  get  bivariate  Poisson  random  variables. 

6 .  Schmeiser's  Bivariate  Gamma  Generator 

Schmeiser  (1979)  developed  a  family  of  algorithms, 
any  member  of  which  can  generate  bivariate  gamma  random 
vectors  having  any  shape  parameters  a^,  and  allowable 
correlation  p.  In  this  section  we  will  discuss  his  general 
procedure . 

If  Z,  and  are  independent  gamma  random  varia¬ 
bles  with  shape  parameters  r,  <5^,  and  6^  >  respectively,  and 
U  is  an  independent  uniform  [0,1]  random  variable,  and  either 
V  =  U  or  V  =  1-U ?  then 

X,  =  F"1^)  +  Z  +  W.  (II-C-5) 

in  1 
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x2  =  Fm1(V)  +  2  +  w2 

where  F  ^  and  F  1  are  the  inverse  distribution  functions 
n  m 

of  gamma  (n,l)  and  gamma  (m,l),  respectively,  distribution 
function.  The  resulting  X^,  x2  are  gamma  random  variables, 
with  shape  parameters 


=  n  +  r  +  5^, 

(II-C-6) 

a2  =  m  +  r  +  5^' 


respectively.  This  result  follows  immediately  from  the 
reproducibility  property  of  the  gamma  distribution  and  from 
noting  that  F  ^(u)  and  F~^(l-u)  are  each  random  variables 
having  CDF  F.  This  scheme  generalizes  the  trivariate 
scheme  but  brings  in  the  inverse  CDF. 

The  correlation  coefficient  is  defined  as 


E[F-1(u) 

n 


P-1 

m 


(v) ]  -  nm  +  r 


X1'X2 


(ala2} 


1/2 


(II-C-7) 


The  remaining  problem  is  to  select  values  of  the  five 
parameters  n,  m,  r,  5^  and  52  to  obtain  the  desired  marginal 
distributions  and  correlation.  The  following  conditions 
must  be  satisfied. 


n  +  r  +  5^  = 
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m  +  r  +  <$2  =  a2 


E[F  1(u)  F  1(v)]  -  (nm)  +  r 
n  m 


(ala2) 


1/2 


=  P 


(II-C-8) 


,  H2,  r,  5lf  6 2  0 


Since  we  are  using  five  variables  to  satisfy  three  equality 
conditions,  finding  a  set  of  parameter  values  corresponds 
to  finding  a  feasible  solution,  rather  than  an  optimal  solu¬ 
tion,  to  a  nonlinear  programming  problem. 

An  efficient  solution  procedure  for  determining 
parameter  values  is  important,  since  substantial  computation 
is  required  to  determine  whether  or  not  conditions  (II-C-8) 
are  satisfied  for  given  parameter  values.  Most  of  the 
computation  is  involved  in  calculating 


p 


E[F"X(u)  F"1(v)] 
n  m 


(c^c^) 


1/2 


nm  -  r 

/ 


since  the  expected  value  must  be  calculated  numerically 
using  any  one  of  the  following  three  integrals: 


/  F~ 1 ( u)  F“1(u)  du  ( II-C-9a) 

0  n  m 


0 


(Fn(x) ) 


xn  exp(-x)  dx/r(n) 


( II-C-9b) 
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(II-C-9C) 


/  Fm1(Fn(-ln  y) ) (-In  y) n  dy/r (n) 

if  p  >  0.  If  p  <  0,  then  replace  Fn^(u)  with  F^Cl-u)  in 
equation  (II-C-9a) ,  replace  Fn(x)  with  1  -  F  (x)  in  equation 
(II-C-9b)  ,  and  replace  Fn(-ln  y)  with  1-F  (-In  y)  in  equation 
(II-C-9c)  . 

Schmeiser  seemed  to  have  best  results  in  terms  of 
a  subjective  trade  off  between  speed  and  accuracy,  using  a 
24  point  Gaussian  method  for  the  integration.  He  selected 
the  parameter  values  from  the  feasible  values  satisfying 
conditions  (II-C-8)  by  making  the  curves  of  regression 
E[XjJX2]  and  E[X2|X^]  behave  as  desired. 


Generating  Procedure 


1. 

Generate  X^  -  gamma 

(n,  1) 

2. 

u  -  F 

n  1 

If  p  <  0,  U  «-  1  -  u 

3. 

Generate  Z  ~  gamma 

(r ,  1) 

4. 

Generate  -  gamma 

(6^1) 

5. 

Generate  ~  gamma 

i — I 

CN 

<o 

6. 

X^  X^  +  z  + 

X2  FH1(U)  +  z  +  w2 
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Based  on  these  procedures  he  developed  a  family  of  algorithms 
which  can  provide  variates  having  any  theoretically  possible 
correlation  p.  Anyway,  for  gamma  marginal  distributions, 
not  all  correlations  are  consistent  with  particular  shape 
parameter  values.  Schmeiser  shows  the  obtainable  correlations 
as  a  function  of  given  =  1  and  5  as  shown  in  Section 

II-A.  Note  that  only  when  =  a ^  is  it  possible  to  obtain 
p  =  1-  Likewise  p  =  -1  is  not  possible  except  in  the  limit 
as  and  a ^  tend  to  infinity.  The  maximum  and  minimum 
possible  correlations,  given  in  Moran  (1969),  occur  when 


X 


2 


and 


respectively,  where  F  (X)  and  F  ^ ( U)  are  the  cumulative 

a  a 


distribution  function  and  inverse  CDF,  respectively,  of  the 
gamma  distribution  with  shape  parameter  a. 
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III.  GENERAL  MIXTURE-TRUNCATION  METHOD 


Denote  by  (Y,Z)  the  bivariate  random  pair,  where  each 
has  identical  marginal  continuous  distribution  F(x) ,  and 
denote  a  general  random  variable  from  this  distribution  by 
X.  The  argument  is  not  specific  to  continuous  random  varia¬ 
bles;  this  aspect  comes  in  only  in  the  computation  of  the 
correlations  and  can  be  developed  in  a  parallel  fashion  for 
discrete  marginal  distributions.  Let 


P  = 


1  - 


1  “  a2  a2 


with  stationary  vector 


1  -  a. 


1  -  a. 


TT  =  77  P  =  ( 


1  —  a^  +  l-  o^7  1  —  a.  +  1  — 


)  = 


( 


and  let  xn  be  an  X  truncated  to  the  left  of  a  fixed  point  x 
1  o 

X  be  an  X  truncated  to  the  right  of  x  ,  so  that 


0 


if  x  <  x 


o 


F  (x) 
x~ 


2 


P[X2ix] 


F  (x)  -  F(xq) 


if  x  >  x 


1  -  F(xo) 


o 


In  addition  we  set  ir^  =  F(xq),  t 


1  -  and  choose 


Y  and  Z  as  follows. 

i)  Choose  Y  from  with  probability  and  then  choose 
Z  from  X^  with  probability  a^,  or  from  X2  with 
probability  1  -  a^. 

ii)  Choose  Y  from  X  w-i-t-h  nrnhahil  ifv  TT  .  Ta/Vi  TT  _  +  TT  =  1, 



and  then  choose  Z  from  X^  with  probability  l-a2,  or 
from  X2  with  probability  a2. 


If  we  choose  (Y,Z)  as  in  the  above  procedure,  then  we  can 
make  the  following  two  theorems. 

A.  MARGINAL  DISTRIBUTIONS 
Theorem  1 . 

The  marginal  distribution  of  (Y,Z)  becomes  F(x)  for  both 
Y  and  Z. 

Proof 

1.  Marginal  distribution  of  Y 

By  definition  Y  is  the  mixture  of  X-^  and  X2  with  proba¬ 
bility  ir 2  respectively.  That  is 


*1  Fx  (x)  +  tt2  Fx  (x) 
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Fy(x)  = 


"l  Flip-  +  ’2  •  0 


F  (x)  -  F(xQ) 

*1  '  1  +  *2  1  -  F ( x  ) 


if  x  <  x 

—  o 


if  x  >  x 


But  since  we  define  tt^  =  F(xq);  =  1  -  ir^  =  1-F(xq),  we 
have 


Fy(x)  = 


F(x  ) 

-7T"  =  F(X) 


F  (x)  -  IT. 


7T,  +  7T 


1  2  IT. 


if  x  <  x 
—  o 


=  F  ( x)  if  x  >  x 


=  F (x)  in  all  cases. 


So  Y  has  the  marginal  distribution  F(x). 

2.  Marginal  distribution  of  Z 
If  Y  is  from  X^,  then 


,  (x)  =  a  F  (x)  +  (l-ai)F  (x) 

‘1  1  X1  X2 


F  (x)  ,  /  \  •  a 

al  F ( X  )  (1  al}  ° 


if  x  <  x 
—  o 


1  +  (1  -  aL) 


F(x)  -  F(xQ) 
1  "  F(xq) 


if  x  >  x 

o 
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If  Y  is  from  X 2,  then 


F  (X)  =  (1-0  F  (X)  +  a-  Fv  (x) 

^2  z  ^  x2 


(1_“2)  fott  +  «2  '  0 


F(x)  -  F(xQ) 

< 1  -  a2>  '  1  +  a2  1  -  F (x  ) 


if 


if 


So, 


Fg  (x)  =  7T 1  Fz  (x)  +  tt ^  F z  (x) 


F (x)  t  ,n  x  F  (x) 
Tt,a,  TTTTT  —  +  ir_  (1  -  cO 


1“1  F (x  )  2  '  2  F (x  ) 

o  o 


F  (x)  -  F(xq) 

7r1(a1  +  (1  -  a±)  1  -  F(X  ) 


if 


if 


F  ( x)  -  F(xq) 
+  TT2((l-a2)  +  a  2  T^fTTT 


and  we  defined  and  tt2  as  follows. 

1  ~  a2 

17 1  1  -  +  1  -  a2 

^2  1  -  +  1  -  (*2 


x  <  x 
—  o 


x  >  x 

o 


x  <  x 

—  o 


x  >  x 

o 
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From  this,  we  know 


"i 

1  -  a  =  -±  (1  -  a,) 

2  n2  1 

If  we  use  this  relationship,  then  F^Cx)  =  F (x)  in  both 
cases.  So  Z  also  has  the  marginal  distribution  F(x).  The 

result  is  a  consequence  of  the  fact  that  tt  is  defined  to  be 

I 

the  stationary  vector  associated  with  P. 

B.  THE  PRODUCT-MOMENT  CORRELATION 
Theorem  2 . 

The  correlation  coefficient  between  Y  and  Z  becomes 

P  =  8  M  , 


where 


-1  <_  8  =  -  (1  -  a 2 )  £  1 


and 


M 


"l"2/ox‘ 


/ 


where 


x 


x  d 


F  (x) 
P(Xo) 
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00 


F  ( x)  -  F  (x  ) 

f  x  d^ - =-7 - r2- 

J  1  -  F  (x  ) 

X  o 

o 


j  x  d 


2  ^  F (x) 


F(xo) 


-  y  • 


2  ^  F(x)  -  F(xQ) 


!  x  d  J  IT 


F (x  )  ”  U2 


X 


/  x2  dF  ( x)  -  E  [X  ] 


Proof 


2  2  2 
=  o ^  tt1  +  o2  tt 2  +  (y1  -  y2)  Tr1  "2 


 cov  [Y ,  Z  ] 


’YZ 


°Y  °Z 


E  [YZ]  -  E [Y]  E  [Z] 
°Y  °Z 


Now 
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E  [YZ] 


where 

Further, 

E[Y]  = 

by  Theorem 


E  [E  [  YZ  |  Y,  Z  e  S]  ] 

w 

E[YZ|Y  e  XirZ  e  X±]P  [Y  e  X± ,  Z  e  X^ 

+  E[YZ|Y  e  X1#Z  £  X2]P[YeX1,Z  e  X2] 

+  E[YZ|YeX2,Z  e  X1]P[Y  £  X2 , Z  £  X±] 

+  E[YZ|Y  £  X2,Z  £  X-jPtY  £  x2,z  £  x2] 

E[X1]  EtX-J  7r1a1  +  ElX^  E[X2]  1^(1  -c^) 

+  E[X2]  E[X1]  ir2(l-a2)  +E[X2]  E[X2]  ^  2  2 
2 

U1  al  *1  +  M1  y2  ^  ~  al^  *1 

+  ^  U  2  ( 1  -  a2)7i2  +  p22  a2Ti2 


M1  =  ElX^,  p2  =  E[X2] 


Es[E[Y|Y  £  S]] 

E[Y|Y  £  XjPtY  £  X±]  +  E[Y|Y  £  X2]P[Y  £  X2] 
E[Xx]  +  E[X2]  tt  2 

yl  ^1  +  p2  "2  = 
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Also  by  Theorem  1, 


2 

a  Y 


If  we  put  together  these  formulae  into 

E[YZ]  ~  E[Y]E[Z] 
PY,Z  ay  a ^ 


E[Y2]  -  (E  [Y]  )  2 


2  2 
ax  °Z 


we  get 


Y ,  Z 


(Ml  -  V2)  Tr1TT2(a1-  ( X  —  ot 2 )  ) 


X 


and  let 


3  =  -  (1  -  a2 ) 

2  2 

M  =  ( M ^  -U2)  Trl7T2//aX 


then 


p  =  3  M  . 

C.  GENERAL  ALGORITHM 

We  give  here  three  algorithms  for  implementing  the 
bivariate  mixture- truncation  method,  which  we  call  the  FXO 
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method,  the  UXO  method,  and  the  TXO  method.  All  of  these 
methods  are  exactly  the  same  except  in  how  the  algorithm 
chooses  xq,  the  truncation  point,  from  the  xq  range  [x^x^]  . 
The  first  procedure,  called  the  FXO,  choose  x  as  a  fixed 

i 

point  of  x^  and  xu  and  uses  the  same  xq  during  the  entire 
routine.  The  second  procedure,  called  the  UXO,  chooses  x 

o 

uniformly  from  [x^,xu]  and  repeats  this  step  in  every  call 
to  the  algorithm.  The  third  procedure,  called  the  TXO,  is 
the  same  as  the  UXO  procedure  except  in  that  it  uses  a  tri¬ 
angular  distribution  instead  of  uniform.  It  is  necessary  to 
fix  these  choices  of  xQ  becuase  in  general  there  is  more  than 
one  xq  which  will  give  a  bivariate  pair  (Y,Z)  with  the  given 
marginal  distribution  and  given  correlation.  The  first 
procedure,  FXO,  is  defective  in  terms  of  their  discontinuity 
of  distribution  while  the  second  and  the  third,  UXO  and  TXO, 
are  satisfactory  in  this  respect.  The  choice  of  the  midpoint 
of  the  interval  [x^,xu]  for  xq  in  FXO  is  based  on  experience 
presented  later  for  various  marginal  distributions.  Note 
that  the  algorithm  described  here  is  inefficient  in  that  it 
generates  the  truncated  variables  X-^  and  X2  by  comparing  ran¬ 
dom  variables  X  to  xq  until  one  which  is  respectively  greater 

than  or  less  than  x  is  found.  More  efficient  methods  can 

o 

be  found  in  special  cases  such  as  the  exponential,  but  the 
present  algorithm  requires  only  a  generation  of  univariate 
random  variables  X  without  regard  to  the  method  used  to  do 
this.  Of  course  initialization  is  required  and  this  is 
specific  to  each  marginal  distribution. 


57 


General  Mixture-Truncation  Method 


1.  (Initialization) 

i)  For  given  marginal  distribution  F(x)  and  correla¬ 
tion  coefficient  p  find  xq  ranges  [x^,xu] 


2.  Define  truncation  point  x 


*  FXO  method 


i)  X, 


2(XJ,  +  Xu) 


*  UXO  method 

i)  Generate  a  uniform  [0,1]  random  variable 
ii)  x  =  x.  +  (x  -  x,  )  *  V. 


I 


u  l. 


*  TXO  method 

i)  Generate  two  uniform  [0,1]  random  variables 

V  V 

ii)  Xq  =  x£  +  xx  +  x2 
where 


x  =  4(x  +  x  ) 

m  2  £  u 


X1  =  (3tm-xt>  *  V1 


X-  =  (x  -  X  )  *  V- 

2  um'  2 


3.  Compute  parameters  value,  ,  it  ,  a^,  a2 


4 .  Choose  type  for  Y 

i)  Generate  a  uniform  [0,1]  random  variable  U 
ii)  If  U  <_  tt^,  go  to  9 
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5.  Y  is  an 

i)  Generate  a  random  variable  X  from  F(x) 

ii)  If  X  >  x  ,  set  Y  X  and  go  to  6 

iii)  Otherwise  return  to  5.i) 

6 .  Choose  type  for  Z 

i)  Set  U  ^  (  (U  -  n1)/(l  -  tt1)  ) 

ii)  If  U  £  1  -  c*2  /  go  to  8 

7.  Z  is  an  X^ 

i)  Generate  a  random  variable  X  from  F(x) 

ii)  If  X  >  x  ,  set  Z  X  and  go  to  11 
o 

iii)  Otherwise  return  to  7.i) 

8.  Z  is  on  X^ 

i)  Generate  a  random  variable  X  from  F(x) 

ii)  If  X  <_  xq/  set  Z  X  and  go  to  11 

iii)  Otherwise  return  to  8.i) 

9.  Y  is  an  X^ 

i)  Generate  a  random  variable  X  from  F(x) 
ii)  If  X  <  xQ,  set  Y  X  and  go  to  10 

iii)  Otherwise  return  to  9 . i ) 

10.  Choose  type  for  Z 
i)  Set  U  «-  U/n1 
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ii)  If  U  <  a^,  go  to  8 
iii)  Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  method  until  a  sufficient 
number  of  random  vectors  are  obtained. 

D.  BIVARIATE  DISTRIBUTION  FUNCTIONS 

From  Theorem  1,  in  Section  III-A,  we  know  that  if  Y 
is  from  X^,  then 


Fz (z |y) 


al  FX  (z)  +  (1 "  “i*  Fx  (z) 


F(z) 


if  z  <  x 


o 


F  ( z )  -  F(xQ) 


if  z  >  x 


o 


and  if  Y  is  from  X  ,  then 


Fz(z|y)  =  (l-a2)  Fx  (z)  +  a2  Fx^  (z) 


if  z  <  x 


o 


F  (  z)  -  F(xQ) 


if  z  >  x 


o 
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By  using  these  we  can  define  the  bivariate  distribution 
function  as  follows. 


F(y,z)  =  P  [Y  <_  y,Z  <_  z] 


I  P[Z  £  z  |  Y 


u]  dP  [Y  <  u] 


where 


P  [Y  <  u]  =  F  (u) 


P  [Z  <_  z  |  Y  =  u] 


,  F(z) 

1  F(xq) 


+  (1  -  a^) 


F  ( z )  -  F(xq) 

-T^rorr 


(1  -  aj 


F(z) 

2 '  F(xo} 


F(z)  -  F(xq) 

(1  -  ot2)  +  a2  l  _  F  (x  ) 


So,  if  we  put  these  together,  integrating  with 


if  u  <  x  ,  z  <  x 

—  o'  —  o 

if  u  <  x  ,  z  >  x_ 

—  o  o 

if  u  >  x  ,  z  <  x 
o  —  o 

if  u  >  x  ,  z  >  x 
o  o 

respect  to 


dP  [Y  <_  u]  ,  we  get  the  final  result: 
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P[Y  <_y,Z  <_z] 


For  example 


-  ^F(z)F(y)  if  Y1V  2-xo  (IH-D-1) 

(l-o,) 

‘  {al  +  rr  [F(z)-F(xq)  ]}F(y)  if  y<xQ,  z>xq  (III-D-2) 

(l-o, ) 

-  (a,  + - — [F(y)-F(x  )  ]  }F(z)  if  y  > x  .  z  <  x^  (III-D-3) 

x  it  2  o  o  —  o 

+  (l-a1)[F(z)-F(xo)]  ^  (III-D-4) 

a2 

+{  (l-a2) +^=-[F(z)-F(xq)  ] }  [F(y)-F(xQ)  ]  if  y>XQ,  z  >  xq 
the  expression  (III-D-4)  is  obtained  as 


F  (y ,  z )  =  /  P[Z  <_  z|Y  =  u]  dF(u) 


V 


j  P[Z  <_  z|Y  =  u  <_  x  ]  dF(u) 

—  00 


+  /  P[Z<_z|Y  =  u>  X  ]  dF(u) 

xo 

(1-0,) 

=  (a,  +  — - - [F  ( z )  -  F(x  )]}  F  ( x  ) 

1  7T  2  O  O 

+  (1-a,)  +  —  [F(z)  -  F  ( x  )  ]  }  [  F  ( y )  -  F(x  _)] 

2  "  2  o  o 

It  is  easily  seen  from  (III-D-2)  that  when  z  », 

P  [Y  <_  y,Z  <_  «]  =  F(y)  ;  from  (III-D-3)  that  as  y  + 
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P  [Y  <_  00 ,  Z  <_  z]  =  F(z)  and  from  (III-D-4)  that  as  y  + 

P[Y  <_  <=°,Z  <_  z  ]  =  F(z)  and  that  z  -»-  °° ,  P  [Y  <_  y,Z  <_  «]  =  F(v)  . 
In  particular  from  (III-D-4)  we  have  that,  as  y  + 


F(y,z)  -*  a1F(xQ)  +  (1  -  a2)  [F(z)  -  F(xq)  ]  +  (l-a2)ir2 
+  a2 [F ( z)  -  F(xq) ] 


=  a1F(xQ)  +  (1  -  a2)  tt2  +  [F  (z)  -  F(x0)] 


=  F (z)  +  (1- 


“2>"2 


-  (1- 


“1>"1 


F(z) 


where  at  the  last  step  we  used  the  facts  that  F(xq)  = 
and  Tr^(l-a^)  =  Tr2(l-a2).  If  F(x)  is  absolutely  continuous 
with  probability  density  function  f(x)  then  the  joint  p.d.f. 
for  the  bivariate  pair  (Y,Z)  is 


f (y#z) 


a,  (1  -  a,  ) 

'  ,2U-  a2)  f(^  f(z) 

i-a 

-  f  (y)  f<*> 

1  -  ai 

-  — — -  f (y)  f(z) 

^2 

a2 

-r  f(y)  f(z) 

*2 


if  y  <  x  ,  z  <  x 
2  —  o  —  o 


lf  y  <  xo. 


if  y  >  xQ,  z  xQ 

if  y  >  x  ,  z  >  x^ 
2  o  o 


(III-D-5) 

( III-D-6 ) 

(III-D-7) 

( III-D-8) 
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Note  that  there  is  a  discontinuity  in  the  density  function 
as  one  crosses  the  boundaries  of  the  four  quadrants  defined 
by  the  lines  y  =  xq;  z  =  xq.  The  density  is  the  same  in  the 
first  and  third  quadrants.  The  multipliers  of  f(y)f(z)/ir2 
are  the  same  in  all  four  quadrants  iff  there  is  independence. 
This  occurs  when  (1-  a^)  =  a 2  so  that  =  a2  and 
f (y fZ)  =  f(y)f(z)  for  the  whole  range  of  y  and  z. 
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IV.  BIVARIATE  UNIFORM  GENERATOR 


The  uniform  random  variable  is  a  continuous  random 
variable  with  probability  density  function  which  is  constant 
over  the  interval  (a,b)  and  zero  otherwise;  the  density 
function 


f  (x)  = 


b-a' 


if  a  <  x  <  b  , 


-  0;  otherwise  ; 


mean  =  E [X]  = 


b+a 
2  ? 


variance  =  VAR[X]  = 


(b  -  a) 
12 


We  note  that  if  U  has  a  uniform  distribution  over  [0,1] 
then  X  =  a  +  (b-a) U  has  a  uniform  distribution  over  [a,b] . 

So  we  will  only  be  concerned  with  algorithms  for  generating 
uniform  [0,1]  distributions. 

Two  algorithms  for  generating  uniform  bivariate  pairs 
are  given  in  Section  II,  Gaver's  and  Lawrance  and  Lewis's, 
and  since  they  have  relatively  smooth  distributed  functions  and 
are  simple  to  generate,  they  are  probably  preferable  to  the 
uniform  bivariate  random  variable  generated  by  the  mixture- 
truncation  method  unless  xq  is  taken  to  be  random  (of  course 
the  correlation  for  the  Gaver  bivariate  uniform  is  difficult 
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to  compute) .  However  the  development  in  this  section  of 
the  uniform  mixture-truncation  method  does  help  to  fix  ideas 
on  the  problem  of  initialization  and  determination  of  the 
possible  range  of  correlations  attainable  in  other  cases. 

To  use  the  mixture-truncation  method  for  generating  bivariate 
random  vector  (Y,Z)  whose  marginal  distribution  is  uniform 
over  [0,1]  identically  and  has  given  correlation  p,  we  should 
use  uniform  [0,1]  distribution  functions  as  F(x).  Then  X^, 
which  is  X  constrained  to  be  on  the  left  side  of  the  trunca¬ 
tion  point  xq,  becomes  uniform  over  [0,xQ]  and  X2 ,  which  is 
constrained  to  be  on  the  right  side  of  x  becomes  uniform 
over  [x  ,1] . 

A.  DETERMINATION  OF  PARAMETERS  IN  THE  UNIFORM  MIXTURE- 
TRUNCATION  METHOD 

Because  X^  is  uniform  (0,xoJ  and  x2  is  uniform  [xq,1], 
E[X1]  =  Ul  = 

VAR[X1]  =  a.2 

e[x2]  =  v2  = 

VAR[X2]  =  a22 


x 

o 

2 


x 


12 


1  +  x 


(l-xo) 

12 


and 


1  -  a. 


=  F(Xo>  =  l-C  +  l-C. 


=  xo  ' 


(IV-A-la) 
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1 


-  a . 


it,  =  1  ~  F(x  ) 

2  O 


1  ”  +  1  -  0*2 


1  -  x  (IV- A- lb) 

o 


If  we  use  these  formulas  in  Theorem  2  of  Section  III,  we  get 


M  =  3  x  (1  -  x  ) 

o  o 


p  =  3  x  (a,  -  x  )  (IV-A-2) 

o  1  o 


I 

Then  from  (IV-A-2),  (IV-A-la)  and  (IV-A-lb) 


a,  =  +  x  (IV-A-3a) 

1  3x  o 


a,  =  1  -  -^(1  -  a.  ) 

2m  TT  2  1 

2  + 

=  1 - - - ^  +  - £ -  (IV- A- 3b) 

1  -  x  3  ( 1  -  x  ) 

o  o 


Furthermore,  and  a ^  are  probabilities  so  they  have  to 
satisfy  the  following  conditions: 


0  <_  <_  1, 

0  _<  a2  ±  1  . 
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By  solving  these  two  inequality  equations,  we  can  get  the 
feasible  range  of  xq  for  given  p  as  follows. 

x.  =  Lower  bound  of  x 

£.  o 

■  1/2  -  /1/4  -(1/3)  p  if  p  >  0 

=  ( IV-A-4) 

/-  p/3  if  p  <  0 

x  =  Upper  bound  of  x 
u  ^  o 

r  1/2  +  v'l/4"-a/3)  P  if  p  >  0 

=  (IV-A-5) 

1  -  /-p/3  if  p  <  0 

After  finding  x„  and  x  ,  choose  x  within  the  range  (x„,x  ] 

by  way  of  either  the  FXO,  UXO,  or  TXO  method.  And  by  formulas 

(IV-A-la) ,  (IV-A-lb) ,  (IV-A-3a)  and  (IV-A-3b) ,  we  can  get 

7T ^ and  c^.  Also  we  can  expect  some  limitations  in 

p  because  of  x  .  From 
o 

M  =  3  x  (1  -  x  ) 

o  o 

-1  <_  3  =  -  (1  -  a 2 )  1 


P  =  6  M 
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We  can  compute  that 


P 


max 


3/4  if  6  =  1,  x,,  <_  xq  =  1/2  <_  xu 


P  • 
mm 


-3/4  if  6  =  -1,  x  <  x  =  1/2  <  x 

a  —  o  —  u 


In  another  way,  note  that  if  one  chooses  xq  =  1/2  and 
c*2  =  =  1,  from  (IV-A-2)  we  attain  the  highest  allowable 

correlation  p  =  3/4,  and  when  =  0,  p  =  -3/4,  the  lowest 
attainable  correlation. 

It  turns  out  that  for  this  choice  of  x  =  1/2  (in  the  FXO 

o 

method) ,  the  bivariate  uniform  distribution  is  also  about 
the  smoothest,  as  can  be  seen  in  a  later  section  (IV-D), 
and  from  the  bivariate  distribution. 

B.  GENERATING  PROCEDURE 

We  developed  here  all  of  the  three  procedures,  the  FXO 
method,  the  UXO  method,  the  TXO  method  for  generating 
bivariate  random  vectors  whose  marginal  distributions  are 
uniform  over  [0,1]  and  correlation  coefficient  is  p.  In 
this  algorithm  we  used  direct  generating  procedures  for  X^ 
and  X2  instead  of  comparing  random  variables  X  to  xq  until 
one  which  is  respectively  greater  than  or  less  than  xQ  is 
found. 
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Uniform  Mixture-Truncation  Method 
1.  (Initialization) 

i)  For  given  -3/4  <  p  <  3/4,  find  x„  and  x 

5c  U 


X 


1/2  -  /1/4  -  1/3  p  if  p  >  0 


p/3 


if  p  <  0 


x 


u 


1/2  +  /1/4  -  1/3  p  if  p  >  0 


1  -  v -p/3 


if  p  <  0 


2.  Define  truncation  point  xq. 


*  FXO  method 


i)  * 


(x,  +  xu)/2 


*  UXO  method 

i)  Generate  a  uniform  (0,1]  random  variable 

ii)  x  =  x  +  (x  -  x  )  *  U, 
o  a  u  z  1 

*  TXO  method 

i)  Generate  two  uniform  [0,1]  random  variables 
Vl'  V2 

111  xo  =  xt  +  X1  +  x2 
where 


X1  =  (xm  '  xi.)  *  V1 
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x0  =  (x  x  )  *  V_ 

2  urn  2 


xm  =  (xu  +  xl)/2 


3.  Compute  tt^,  ■n^,  and  a 2  as 


"l  =  F(xo>  =  xo 


"2  ’  1  ‘  "l 


^  X  +  X 

3  o  o 


1  -  -±(1  -  a  ) 
11 2  1 


4.  Choose  type  for  Y 

i)  Generate  a  uniform  [0,1]  random  variable  U 


ii)  If  U  <_  ir^  go  to  9 


5.  Y  is  on  X2 

i)  Generate  a  uniform  [0,1]  random  variable  W1 


ii)  Y  x  +  (1.0  -  x  )  *  W. 

o  o  1 


6 .  Choose  type  for  Z 

i)  Set  U  ^  (  (U  -  tt1)/(1  -  ) 


ii)  If  U  <  1  -  a.,  go  to  8 
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7. 


Z  is  an  X2 


i) 

Generate 

a  uniform 

ii) 

Z  x  + 

o 

(1.0  -  xc 

iii) 

Go  to  11 

Z 

is  an  X^ 

i) 

Generate 

a  uniform 

ii) 

Z  <-  x  * 
o 

w2 

iii) 

Go  to  11 

Y 

is  an  X^ 

i) 

Generate 

a  uniform 

ii) 

Y  «-  x  * 
o 

W1 

[0,1]  random  variable  W 2 


[0,1]  random  variable  W2 


[0,1]  random  variable 


10.  Choose  type  for  Z 

i)  Set  U  «-  U/-n1 

ii)  If  U  <  go  to  8 

iii)  Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4,  for  the  FXO  method,  or  go 
to  2  for  the  UXO  method  and  the  TXO  method  until  a 
sufficient  number  of  random  vectors  are  obtained. 


The  programs  are  listed  in  the  appendix  and  the  scatter 
)lot  and  bivariate  histogram  of  resulting  random  vectors  are 
>hown  in  Section  IV-D. 
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c.  REGRESSION  OF  Z  ON  Y  FOR  GIVEN  p 


Schmeiser  (1979)  in  developing  a  bivariate  gamma  method 

fixed  the  free  parameter  by  specifying  the  regression  of  Z 

and  Y  and  vice-versa.  We  do  this  now  for  the  uniform  case 

with  fixed  x  and  uniformly  distributed  x  .  First  in  fixed 
o  o 

x  case  we  have,  for  y  <  x 
o  —  o 


E  [Z  I  Y  =  y,y<xo] 


E[X]_]  +  (1  -  ax)  E[X2] 


and  for  y  >  x 

o 


E[Z|Y  =  y,y  >  Xq]  =  (1  -  a2)  E[X1]  +  a2  E[X2] 


These  are  the  step  function  and  clearly  not  linear  in  y. 

If  x  has  uniform  distribution  over  [x.,x  ],  we  have  the 

O  X)  u 

following  results.  These  results  are  dependent  on  Y  values. 
If  y  <_  x^,  then  we  have 


E  [Z  |  Y  =  y  ]  =  /  E  [Z  |  Y  =  y  ,X  =  XQ,y  <_  xQ]  f  (xq)  dxQ 


x 


£ 
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E[Z|Y  =  y] 


- - - (~  X  -  \  X  -  £  In  — ) 

xu "  xi  2  u  2  *  6  x,i 


If  x»  <  y  <  x  ,  then  we  have 

x,  —  —  U 


E  [Z  |  Y  =  y]  =  J  E  [  Z  |  Y  =  y ,  X  =  xo  ,  y>xQ]f(xo)  dxQ 


X* 


X 


u 


+  j  E[z|Y  =  y,  X  =  x  ,  y<x  ]  f(x  )  dx 

y 


sv  4  + 


2  •  etT^T1  f(xo>  dxo 
x.  o 


X 


w  1 

+  j  (j 

y 


f(x  )  dx 
6x  o  o 


1  ,1  1  p  ,  1-y  xu  , 

-  (tT  X  -  TT  X  -7  In  — *-  ■= - ) 

x  -  x^  2  u  2  1  6  y  1-x 


If  y  >  x  /  then  we 


x 

u 

E  [  Z  |  Y  =  y  ]  =  j  E  [Z  |  Y  =  y ,  X=xQ,  y  >  xQ]  f(xQ)  dxc 

XA 


u  1 

/  <7  + 


I  -r-f-, - - - r-)  f(x  )  dx 

2  6(1-  xQ)  o'  o 


1  ,1  1  p  ,  1  xu, 

(n-  X  -r  X,  -  7  In  i— — ) 


x-x2u2£6  1-x 

U  A  £ 
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From  the  results  we  can  see  that  if  we  use  uniform  distri- 

bution  for  xq  we  can  get  smoother  regression  functions  if 

x-o  1  Y  £  xu  although  we  still  have  step  functions  for  the 

y  <  xfl  and  y  >  x  cases. 

1  —  H  1  —  u 

D.  SIMULATION  RESULTS 

We  will  show  here  the  scatter  plots  and  the  bivariate 
histograms  of  the  mixture-truncation  bivariate  uniform  (0,1) 
random  vectors  with  correlation  p  =  0.3  and  p  =  -0.3  by  the 
FXO,  UXO  and  TXO  methods  in  Figures  (IV-a),  (IV-b)  and 
(IV-c) ,  respectively. 

In  addition,  we  will  show  the  scatter  plot  and  bivariate 
histogram  of  sample  size  2000  from  Gaver's  transformation  and 
Lawrance-Lewis  method  in  Figures  (IV-e)  and  (IV-d) , 
respectively.  The  Gaver  method  does  not  have  known  correla¬ 
tion  but  the  value  of  p  to  give  p  in  the  exponential  case 
is  used. 

As  we  mentioned  earlier,  the  FXO  method  has  some  discon¬ 
tinuity  at  the  truncation  point.  But  the  UXO  and  TXO  methods 
generate  relatively  smooth  continuous  distributions.  In  this 
respect,  we  say  that  the  FXO  is  defective  and  the  UXO  and 
the  TXO  are  relatively  satisfactory. 

The  computed  correlation  from  subroutine  BIVHST  (Bivariate 
histogram)  is  a  little  different  from  given  correlation.  But 
this  can  be  assumed  as  a  sampling  error.  To  check  this  error, 
we  simulated  10  times  with  given  correlation  p  =  0.1.  From 
this  simulation  we  get  the  following  results: 
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p  =  mean  of  computed  correlation  =  0.105 


VAR[p]  =  variance  of  computed  correlation  = 


o[p]  =  standard  deviation  of  mean  =  0.0064 


0.0004 
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Figure  IV-al.  Scatter  plot  for  sample  of  size  2000  from 

mixture-truncation  bivariate  uniform  with 
p  =  0.3.  Here  xQ  is  fixed  at  the  midpoint 
between  the  lower  and  upper  bounds. 
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Figure  IV-a2. 


Bivariate  historgram  for  sample  of 
from  mixture-truncation  bivariate 
with  p  =  0.3.  Here  x  is  fixed  at 
midpoint  between  the  ?ower  and  the 
bounds . 
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Figure  IV-a3.  Scatter  plot  for  sample  of  size  2000  from 

mixture- truncation  bivariate  uniform  with 
p  =  0.3.  Here  x  is  fixed  at  the  midpoint 
between  the  lower  and  the  upper  bounds. 
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Figure  IV-a4 . 


Bivariate  histogram  for  sample  of  size  2000 
from  mixture-truncation  bivariate  uniform 
with  p  =  -0.3.  Here  x  is  fixed  at  the 
midpoint  between  the  l8wer  and  the  upper 
bounds . 
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Figure  IV-bl.  Scatter  plot  for  sample  of  size  2000  from 

mixture-truncation  bivariate  uniform  with 
p  =  0.3.  Here  x  is  taken  to  be  uniformly 
distributed  between  the  lower  and  the 
upper  bounds . 
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Figure  IV-b2 . 


Bivariate  histogram  for  sample  of 
from  mixture-truncation  bivariate 
with  p  =  0.3.  Here  x  is  taken  to 
uniformly  distributed  between  the 
the  upper  bounds . 
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Figure  IV-b3.  Scatter  plot  for  sample  of  size  2000  from 

mixture-truncation  bivariate  uniform  with 
p  =  -0.3.  Here  xq  is  taken  to  be  uniformly 
distributed  between  the  lower  and  the  upper 
bounds 
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Figure  IV-b4. 


Bivariate  histogram  for  sample  of  size  2000 
from  mixture-truncation  bivariate  uniform 
with  p  =  -0.3.  Here  x  is  taken  to  be 
uniformly  distributed  Between  the  lower  and 
the  upper  bounds. 
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Figure  IV-cl.  Scatter  plot  for  sample  of  size  2000  from 

mixture-truncation  bivariate  uniform  with 
p  =  0.3.  Here  xQ  has  triangular  distribution 
between  the  lower  and  the  upper  bounds . 
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Figure  IV-c2. 


Bivariate  histogram  for  sample  of  size  2000 
from  mixture-truncation  bivariate  uniform 
with  p  =  0.3.  Here  x  has  triangular 
distribution  between  the  lower  and  the 
upper  permissible  bounds. 
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Figure  IV-c3.  Scatter  plot  for  sample  of  size  2000  from 
mixture- truncation  bivariate  uniform  with 
p  =  -0.3.  Here  x  has  triangular  distribu¬ 
tion  between  the  Sower  and  the  upper  bounds  . 
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Figure  IV-c4 . 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
uniform  with  p  =  -0.3.  Here  x  has 
triangular  distribution  between  the  lower 
and  the  upper  bounds . 
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Figure  IV-dl . 


Scatter  plot  for  sample  of  size  2000  from 
Lawrance-Lewis ' s  method  with  p  =0.3  and 
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Figure  IV-d2. 


Bivariate  histogram  for  sample  of 
from  Lawrance-Lewis ' s  method  with 
and 
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Figure  IV-el. 


Scatter  plot  for  sample  of  size  2000 
from  Gaver ' s  transformation  with  p  =  0.6. 
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Figure  IV-e2.  Bivariate  histogram  for  sample  of  size 

2000  from  Gaver ’ s  transformation  with 

p  =  0.6. 
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V.  BIVARIATE  EXPONENTIAL  GENERATOR 


Another  probability  distribution  of  major  interest  in 
simulations  is  the  exponential.  The  cumulative  distribution 
function  and  probability  density  function  for  the  exponential 
are,  respectively. 


F(x)  =  l-e^X  x>0 


and 


f  (x)  =  A  e  Ax  x  >_  0 

The  expected  value  of  the  exponential  distribution  is 

E [X]  =  1/A 


and  the  variance  is 


VAR [X]  =  1/A2  . 

The  problem  of  generating  exponential  deviates  reduces  to 
one  of  generating  "unit"  exponentials,  i.e.,  those  with 
A  =  1.0,  and  then  multiplying  the  result  by  whichever  1/A 
is  necessary  to  give  the  desired  distribution.  That  is,  if 
the  random  variable  E  has  the  exponential  (A  =  1)  distribution 
then, 
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X 


las  the  exponential  (A  =  n)  distribution.  In  this  section, 
tfe  considered  only  unit  exponentials  as  marginal  distributions 
for  bivariate  pairs. 

DETERMINATION  OF  PARAMETERS  IN  THE  EXPONENTIAL  MIXTURE- 

TRUNCATION  METHOD 

Because  X-^  is  an  exponential  (A  =  1)  truncated  to  the 
left  of  xQ  and  X2  is  also  exponential  (A  =  1)  truncated  to 
the  right  of  xq/  we  have 


E[X1] 


0 


1 


VARfX-J 


2 


0 


E[X2] 


oo 


F(x)  -  F(xo) 


X 


o 


1  +  x 


o 
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VAR [X2 ] 


x 


00  0  F  (x)  -  F  (x  ) 

/  1-F(x°) 


-  y . 


=  1; 


ind  from  definition. 


1  -  a. 


1  1  -  +  1  -  0.2 


-  P<Xo> 


-x 


=  1  -  e 


(V-A-la) 


1  -  a . 


2  1  -  +  1  -  0.2 


=  1  -  F(xq) 


-x 


=  e 


(V-A-lb) 


If  we  use  these  formulas  in  Theorem  2  in  Section  III,  we  get 


M  = 


2  2 
—  x 

7T^  O 


(V-A-2) 


3  =  — (a,  -  it,) 

7T  2  1  1 


(V-A-3) 


and 
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al  -  "i  2 

p  —  -  x 

TT  ,  O 


( V-A-4 ) 


For  given  correlation  coefficient  p,  we  can  compute  as 
a  function  of  xQ  from  the  formula  (V-A-4)  as 


P  IT 


\  +  ^1 


-x 


=  (1  +  -£t)  (1  -  e  °) 

x 


and  we  know  that  a0  =  1 - (1  -  a, )  from  the  formulas  (V-A-la) 

2  it  2  1 

and  (V-A-lb) .  From  this, 


a2  =  1  "  ^“‘l* 


-x  x  -x  n 

=  e  °  +  e  °<1  -  e  °)2  -A 

X 


These  and  a ^  are  probabilities,  so  they  have  to  be  greater 
than  or  equal  to  zero  and  less  than  or  equal  to  one; 


0  < 


<  1 


0  <  e  °  +  eX°(l 


-x 


Oj  2 


< 


1 


(V-A-5a) 

(V-A-5b) 


From  these  two  inequality  equations,  (V-A-5a)  and  (V-A-5b) , 
we  can  find  the  xq  ranges  for  given  correlation  coefficient 
p.  To  solve  these  equations,  we  can  divide  into  two  cases. 
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one  for  positive  correlation  case  and  the  other  for  negative 
correlation  case.  If  correlation  coefficient  p  is  positive, 
then  both  equations  are  always  positive.  Thus  we  only  need 
to  find  the  xQ  ranges  which  makes  (V-A-5a)  and  (V-A-5b)  less 
than  1.  From  the  equation  (V-A-5a) ,  for  the  ct.  case, 

-x 

ax  =  (1  -  e  °)  (1  +  •)  <  1 

X 

o 

becomes 

x  n 

p  (e  -  1)  <.  x 

and  let 

~  x 

Z  /  O  T  X 

Y1  =  XQ  /  y2  =  p  (e  -  i)  • 


Because  of  the  first  derivatives  of  y^  and  Y2  are  always 
positive,  we  know  that  these  two  functions  are  monotone 
increasing  functions.  Thus  we  can  find  xQ  ranges  which 
satisfy  y^  >_  Y2  by  the  Newton  Raphson  method.  When  using 
the  Newton  Raphson  method,  let  Y  ~  Y]_  ~  Y2  =  0  and  find  an 
approximate  solution,  by  approximating  exponential  series, 
which  we  can  use  as  a  starting  point.  That  is. 


x 


x 


y  -  y1-y2 


xo  '  0(1  +  xo  +  TT 


3! 


-  1) 


x  + 
o 


0 
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Then 


x 

o 


U-t>  1  (1  "  p  ~  TZ 


p2)1/2 


1/3  p 


Starting  with  this  approximate  value  in  the  Newton  Raphson 
method,  we  can  find  xQ  range,  say  [x„  ,  x^  ] ,  which  satisfies 

0  <_  <_  1.  And  for  the  case, 

-x  -x  -x  ~ 

a 2  =  e  °  +  e  °(l-e  °)2  1  1 

xo 

becomes 


x  ~ 

p  (e  -  1)  <  x 

—  o 

This  result  is  exactly  the  same  as  the  case,  that  is,  at 

the  same  range  a-^  and  satisfy  constraints  <_  1  and 

a_  <  1.  Thus  we  can  use  x.  and  x  as  the  lower  bound  of 

2  -  *i  U1 

x  .x  ,  and  the  upper  bound  of  x  ,x  .  If  correlation  coeffi- 
o  u'  ^  o  u 

cient  p  is  negative,  then  the  equations  (V-A-5a)  and  (V-A-5b) 
are  always  less  than  1.  Therefore  we  need  to  consider  only 
one  constraint  which  makes  _>  0,  a2  >.  0.  From  the 
equation  (V-A-5a) ,  solve  the  inequality  equation 

0  <_  «i  =  (1  -  e  X°)  (1  +  -2^) 
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-x 

Since  1  -  e  °  is  always  positive,  we  see  that  to  satisfy 
the  inequality  1  +  — ^  should  be  positive,  i.e., 


P 


2 


x 


o 


> 


-1 


equivalently,  we  have 


x 


o 


> 


In  the  (*2  case,  from  equation  (V-A-5b)  , 

-x  x  -x  n 

0  <  a2  =  e  °  +  e  °(1  -  e  °)  2 

xo 

or,  equivalently,  we  have 

o  x  „ 

xQ  _>  -p  (e  -  1) 

As  in  the  positive  correlation  case,  we  can  find  a  starting 
point  by  approximation  to  solve  this  equation  by  Newton 
Raphson.  The  result  comes  out  as 


Xg  =  (/^7  +  p)/(-  j) 

with  this  starting  point  we  find  another  bound  of  xq  which 
satisfies  0  <  a-.  This  becomes  the  upper  bound  of  x  ,  x  , 
and  from  the  case,  we  have  a  constraint  xq  >_  f-p  which 
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becomes  the  lower  bound  of  x  ,  i.e.,  x  =  /-p.  The  lowest 

and  highest  correlations  available  for  bivariate  exponential 

pairs  in  mixture-truncation  method  are  approximately  -0.480  and 

0.647  respectively.  By  comparison  note  that  the  most  negative 

correlation  available  for  bivariate  exponential  pairs  with 

2 

identical  fixed  marginals  is  1-^-  [Moran  (1967)].  Gaver's 

(1972)  negatively  correlated  pair  has  correlations  in  the 

range  [-0.5,0].  The  table  (V-l)  shows  the  lower  and  upper 

bound  of  x  in  the  mixture-truncation  method  with  identical 
o 

marginal  exponential  and  given  correlation. 


Table  V-l:  The  lower  bound  and  upper  bound  of  x 
in  the  mixture- truncation  method  witn 
identical  exponential  marginal  distributions 


xq  range  for  each  correlation 


p 

XL 

X 

u 

P 

XL 

X 

u 

0.1 

0.106 

5.832 

-0.1 

0.317 

1.984 

CM 

• 

o 

0.225 

4.723 

CM 

• 

O 

1 

0.448 

1.439 

0.3 

0.362 

3.990 

1 

o 

• 

u> 

0.548 

1.103 

0.4 

0.527 

3.395 

• 

o 

1 

0.633 

0.855 

0.5 

0.741 

2.842 

-0.45 

0.671 

0.751 

0.6 

1.  082 

2.223 

B.  GENERATING  PROCEDURE 

The  generating  procedure  of  bivariate  exponential  random 
vector  is  almost  the  same  as  .in  the  uniform  case.  As  in 
the  uniform  case,  we  develop  three  procedures  which  are  the 
FXO  method,  the  UXO  method  and  the  TXO  method.  As  mentioned 
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in  Section  III,  all  of  these  methods  are  exactly  the  same 

except  in  how  xq  is  chosen  from  the  xQ  range  [x^,xu].  The 

FXO  method  chooses  xQ  as  a  fixed  point  between  x,  and  x^. 

This  procedure  makes  some  discontinuity  at  the  truncation 

point.  In  this  respect,  this  method  is  defective. 

The  UXO  and  TXO  methods  are  the  same  except  we  use 

different  distributions  for  x  .  In  the  UXO  method  we  use 

o 

the  uniform  distribution  and  the  triangular  distribution  for 
the  TXO  method.  These  methods  have  more  smooth  distribution 
than  the  FXO  method. 


Exponential  Mixture-Truncation  Method 
1.  (Initialization) 

i)  For  given  -0.48  <  p  <  0.64,  find  x^  and 


2.  Define  truncation  point  xq 

*  FXO  method 

i)  x„  =  i(x„  +  x  ) 
o  2  £  u 

*  UXO  method 

i)  Generate  a  uniform  [0,1]  random  variable  U^ 


ii)  xo  =  x,  +  (xu  -  xt)  *  U1 
*  TXO  method 

i)  Generate  two  uniform  [0,1]  random  variables 
V  V2 

ii)  X0  =  X£  +  X 1  +  X2 


where 
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xra  =  (x*  +  V/2 


x. 


=  ( x  - 


m 


V 


V, 


=  (X  - 


u 


xm> 


V, 


3.  Compute  parameter  values 

-x 

tt  ^  =  F  (xq)  =  1  -  e  ° 

"2  “  1  *  ”l 

“l  =  "l(1  +  ~£1) 

xo 

a0  =  1 - -{1  -  a,) 

2  ^2  1 

4.  Choose  type  for  Y 

i)  Generate  a  uniform  [0,1]  random  variable  U 

ii)  If  U  £  go  to  9 

5.  Y  is  an 

i)  Generate  an  exponential  random  variable 
ii)  If  E^  >  xq,  set  Y  +■  E^  and  go  to  6 

iii)  Otherwise,  return  to  5.i) 

6.  Choose  type  for  Z 

i)  Set  U  ^  ( ( U  -  tt1)/(1  -  tt1)  ) 
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ii)  If  U  <  1  -  a^,  go  to  8 

7.  Z  is  an  Xj 

i)  Generate  an  exponential  random  variable 
ii)  If  E2  >  xq/  set  Z  E2  and  go  to  11 
iii)  Otherwise  return  to  7.i) 

8.  Z  is  an 


i) 

Generate  an  exponential  random  variable  E2 

ii) 

If  E2  <_  xq/  set  Z  *■  E2  and  go  to  11 

iii) 

Otherwise  return  to  8.i) 

9 .  Y  is  an  X:^ 


i) 

Generate  an  exponential  random  variable  E^ 

ii) 

If  E^  <_  xQ,  set  Y  E^  and  go  to  10 

iii) 

Otherwise  return  to  9  .  i ) 

10.  Choose  type  for  Z 


i) 

Set  U  U/t^ 

ii) 

If  go  to  8 

iii) 

Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  methods  until  a  sufficient 
number  of  random  vectors  are  obtained. 

For  the  exponential  case  it  is  possible  to  give  a  more  effi 
cient  algorithm  in  which  X^  and  X2  are  generated  exactly. 
The  algorithm  is  as  follows. 
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Efficient  Exponential  Mixture-Truncation  Method 


1.  (Initialization) 

i)  For  given  -0.48  <  p  <  0.64,  find  xp  and  xu 

2.  Define  the  truncation  point  xq 

*  FXO  method 

1/  ^  X 

X  =  T  X.  +  X  ) 

o  2  1  u 

*  UXO  method 

i)  Generate  a  uniform  [0,1]  random  variable 


ii)  x 


o 


*  TXO  method 

i)  Generate  two  uniform  [0,1]  random  variables 


ii)  x 


o 


where 


x 


m 


x 


1 


x 


2 


3.  Compute  parameter  values 


1  -  e 


o 
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1 


n 


2 


al  =  *1^  +  ~^~2^ 

xo 

“2  "  1  -  ^(1  -  al> 

4.  Choose  type  for  Y 

i)  Generate  a  uniform  [0,1]  random  variable  U 
ii)  If  U  <  tt^,  go  to  9 

5.  Y  is  an 

i)  Generate  an  exponential  random  variable 

ii)  Set  Y  x  +  E, 
o  1 

6.  Choose  type  for  Z 

i)  Set  U  ^  ( (U  -  tt1)/(1  -  ir1)  ) 

ii)  If  U  <  1  -  o.^,  go  to  8 

7.  Z  is  an  X^ 

i)  Generate  an  exponential  random  variable  E ^ 

ii)  Set  Z  -ex  +  E-  and  go  to  11 
o  2  3 

8 .  Z  is  an  X^ 

i)  Generate  a  uniform  [0,1]  random  variable 
ii)  Set  Z  -  ln(1.0  -  V^*^)  and  go  to  11 

9.  Y  is  an  Xx 

i)  Generate  a  uniform  [0,1]  random  variable 
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ii)  Set  Z  -  ln(1.0  - 

10.  Choose  type  for  Z 
i)  Set  U  4-  U/*1 

ii)  If  U  <  a^,  go  to  8 

iii)  Otherwise,  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  methods  until  a  sufficient 
number  of  random  vectors  are  obtained. 

Note  that  to  compute  and  in  step  1  of  both  algorithms 
we  use  subroutine  BOUND  which  is  listed  in  the  appendix  and 
the  program  with  effective  algorithm  also.  The  scatter  plot 
and  bivariate  histogram  of  resulting  random  vectors  are  shown 
in  Section  V-D. 

C.  REGRESSION  OF  Z  ON  Y  FOR  GIVEN  p 

Schmeiser  (1979)  has  used  the  regression  of  Z  on  Y = y 
to  fix  the  parameters  in  his  bivariate  gamma  distribution. 
Consequently  we  investigate  this  for  the  mixture-truncation 
method  case.  The  regression  is  different  depending  on 
whether  Y  <  x^  or  Y  >  x  .  We  consider  two  cases  here,  one 
for  fixed  xQ  and  the  other  for  xq  having  uniform  distribution. 
For  fixed  xq,  we  have 


E[Z|Y  =  y,y  <_xq] 


a1E[x1]  +  (l-a1)E[x2] 


1  +  x 


o 
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Substituting  the  value  for 


a 


1 


then  we  have 

E[Z|Y  =  y,y  <  xq]  =  1  +  xq  -  xq  ( 1  +  -^) 


And  if  y  >  xq/  then 

E[Z|Y  =  y,y  >  xq]  =  ( 1  -  a2 )  E  ]  +  c^E  [x2  ] 


Substituting  the  value  for 


a 


2 


then  we  have 


E  [Z  |  Y  =  y  ,y  >  xq] 


1  +  —  di_ 

TT ,  X 
2  O 
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Thus  the  regression  is  constant  over  (0,xq)  and  changes  for 
y  >_  xq .  This  is  not  surprising  in  light  of  the  joint  dis¬ 
tribution  given  in  Section  IV-D.  For  uniformly  distributed 
xq/  the  computation  is  different  for  different  ranges  of  Y. 
If  y  <_  x  ,  then  we  have 


x 

u 

E  [Z  |  Y  =  y  ]  =  /  E  [Z  |  Y  =  y  /  X  =  xq  ,  y<XQ]  f(xQ)  dxc 


If 


X 


u 


X 


I  d - E2)  e  °  dx0 


-X  -x. 


x  -x 

--U  (  u  e  ° 

-  e  +  e  -  p  j  - -  dx 


x 


2  o 


£ 


-XJl  -XU  ,  ,  1  “XU 

e  ■  e  +  p  (“  e 

u 


1  "X£ 
—  e 

X£ 


X 

+  In  —  -  x  +  x. ) 
x„  u  £ 


x„  <  y  <  x 
£  —  J  —  u 


then  we  have 
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E  [Z  |  Y  =  y  ]  =  J  E  [Z  |  Y  =  y ,  X=xq,  y>_XQ]f(xo) 

*L 


dx 


x 


u 


+  /  E  [Z  |  Y  =  y ,  X=xq,  y<_xQ]f(xo)  dxc 

y 


=  2e  y  +  2py  -  p (xp  +  x  )  +  p ( 


-x 


+  x  )  +  p  (-i-  e  u  -  y 


i,  u  "  'x  y 

u  J 


x 

+  p  In  — 

y 


If  y  >  x  ,  then  we  have 
1  u 


X 

u 

E  [Z  |  Y  =  y  ]  =  I  E  [Z  |  Y  =  y  ,  X  =  Xq,  y>XQ]f(xo)  dXQ 

x  „ 


X 


U  7T,  -X 

J  (1  +  —  — )  e  °  dx 

J  TT-  X  O 

2  o 


""X  —  x 

e  £-e  u  +  p  (y  ~  x  ) 


By  making  xq  uniformly  distributed  over  the  available  range 
of  xq  for  given  correlation,  we  can  get  smoother  behavior 
for  the  regression  function. 

D.  SIMULATION  RESULTS 

We  will  show  here  the  scatter  plots  and  bivariate 
historgram  of  the  mixture-truncation  bivariate  exponential 
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random  vectors  with  correlation  p  =0.3  and  p  =  -0.3  by 
the  FXO,  UXO  and  TXO  methods  in  Figures  (V-a) ,  (V-b)  and  (IV-c) , 

respectively.  Also  we  show  the  results  from  Gaver's  bivariate 
exponential  and  Marshall  and  Olkin' s  in  Figures  (V-d)  and 
(V~e) ,  respectively. 

In  the  mixture-truncation  bivariate  exponential  distribu¬ 
tion,  the  generated  random  vectors  by  the  FXO  method  also 
has  discontinuity  at  truncation  point  but  the  other  cases 
have  relatively  smooth  distributions. 

The  computed  correlation  printed  in  the  left  side  of  the 
bivariate  histogram  is  a  little  different  from  the  given 
correlation.  But  this  difference  can  be  assumed  as  a 
sampling  error. 

To  check  this  error  we  simulated  10  times  with  given 
correlation  p=  0.1.  From  these  simulations  we  get: 

p  =  mean  of  computed  correlation  =  0.092 

VAR[p]  =  variance  of  computed  correlation  =  0.0008 

a(p)  =  standard  deviation  of  mean  =  0.009. 
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X  X 


Figure  V-al.  Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  0.3.  Here  xq  is  fixed  at  midpoint 
between  the  lower  and°the  upper  bounds  . 
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Figure  V-a2 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
exponential  with  p  =  0,3.  Here  x  is 
fixed  at  the  midpoint  between  the°lower 
and  the  upper  bounds . 


112 


Figure  V-a3.  Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  -0.3.  Here  x  is  fixed  at  the 
midpoint  between  the  lower  and  the  upper 
bounds. 
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Figure  V-a4 . 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
exponential  with  p  =  -0.3.  Here  x  is 
fixed  at  the  midpoint  between  the  Sower 
and  the  upper  bounds . 
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Figure  V-bl.  Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  0.3.  Here  x  is  taken  to  be 
uniformly  distributeabe tween  the  lower 
and  the  upper  bounds  . 
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Figure  V-b2 .  Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
exponential  with  p  =  0.3.  Here  x  is  taken 
to  be  uniformly  distributed  between  the 
lower  and  the  upper  bounds. 


116 


*C.  rM€P\M|C(iS  0  li  23  31  39  4?  33  63  71  79  67  93  103  111  119  |*T 


Figure  V-b3.  Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  -0.3.  Here  x  is  taken  to  be 
uniformly  distributed  between  the  lower 
and  the  upper  bounds  . 
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Figure  V-b4. 


Bivariate  histogram  for  sample  of  size  2000 
from  mixture-truncation  bivariate  exponential 
with  p  =  -0.3.  Here  x  is  taken  to  be 
uniformly  distributed  between  the  lower  and 
the  upper  bounds  . 
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Figure  V-cl.  Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  0.3.  Here  x  has  triangular 
distribution  between  ?he  lower  and  the 
upper  bounds. 
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Figure  V-c2. 


Bivariate  histogram  for  sample  of  size  2000 
from  mixture-truncation  bivariate 
exponential  with  p  =  0.3.  Here  x  has 
triangular  distribution  between  tfte  lower 
and  the  upper  bounds . 
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Figure  V-c3.  Scatter  plot  for  sample  of  size  2000  from 
mixture-truncation  bivariate  exponential 
with  p  =  -0.3.  Here  x  has  triangular 
distribution  between  tRe  lower  and  the 
upper  bounds. 
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Figure  V-c4. 


Bivariate  historgram  for  sample  of  size 
2000  from  mixture-truncation  bi¬ 
variate  with  p  =  -0.3.  Here  x  has 
triangular  distribution  between  the  lower 
and  the  upper  bounds. 
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Figure  V-dl.  Scatter  plot  for  sample  of  size  2000  from 
Gaver ' s  bivariate  exponential  with 
P  =  -0.3. 
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Figure  V-d2 .  Bivariate  histogram  for  sample  of  size 
2000  from  Gaver's  bivariate  exponential 
with  p  =  -0.3. 
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Figure  V-el.  Scatter  plot  for  sample  of  size  2000  from 
Marshall  and  Olkin's  bivariate  exponential 
with  p  =  0.3. 
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Figure  V-e2.  Bivariate  histogram  for  sample  of  size  2000 
from  Marshall  and  Olkin ' s  bivariate 
exponential  with  p  =  0.3. 
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VI.  BIVARIATE  GAMMA  GENERATOR 


The  gamma  distribution  with  shape  parameter  r  and 
scale  parameter  X  has  the  density  function 


f  (x)  = 


r  (r) 


r-1  -Xx 
x  e 


if  x  >  0 


if  x  <  0 


(VI-0-1) 


and  cumulative  distribution  function 


F(y)  =  j 


y  xr 


r  (r) 


r-1  -Xx  , 
x  e  dx  . 


(VI-0-2a) 


In  particular,  if  the  shape  parameter  r  is  a  positive  integer 
then 


F(y) 


1 


r-1 


l 

j=0 


e  Ay(X y) 

j! 


where  r (r)  is  Euler's  gamma  function 


(VI-0-2b) 


T (r)  =  /  xr  ^  e  x  dx  . 

0 

If  the  random  variable  X  has  Gamma  (r,X)  distribution  then 
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E  [x) 


r 

X 


VAR [x]  =  -% 

when  r  =  1,  X  has  the  exponential  distribution  while  X, 

suitably  scaled,  has  an  asymptotically  normal  distribution 

as  r  ->  ».  We  note  that  if  X  has  a  Gamma  (r,l)  distribution 

x 

then  —  has  a  Gamma  (r,X)  distribution.  So  we  may  set  =  1 
as  far  as  the  generating  algorithm  is  concerned.  The  output 
from  the  generator  may  then  be  appropriately  scaled.  Two 
algorithms  for  generating  bivariate  gamma  random  vectors  are 
given  in  Section  II.  Since  they  have  complicated  computa¬ 
tions  to  define  the  parameter  values  and  use  inverse  trans¬ 
formation  functions,  the  mixture-truncation  method  is 
probably  preferable  to  those  methods.  To  generate  bivariate 
random  vectors  with  identical  Gamma  (r,l)  marginal  distribu¬ 
tions  and  having  given  correlation  p  by  mixture-truncation 
method,  we  only  need  univariate  gamme  (r,l)  random  variate 
generator.  Under  the  assumption  that  the  univariate  gamma 
generator  is  available  we  developed  this  procedure.  In  fact 
we  used  Lewis  and  Robinson's  gamma  generator  as  a  univariate 
gamma  generator. 

A.  DETERMINATION  OF  PARAMETERS  IN  THE  GAMMA 
MIXTURE -TRUNCATION  METHOD 

If  the  marginal  distribution  is  gamma  with  positive 
integer  shape  parameter  r  and  scale  parameter  X  =  1,  then 
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we  know  that  F(x)  becomes  cumulative  distribution  function 

of  gamma  (r,l).  Then  the  is  gamma  (r,l)  truncated  to  the 

left  of  xq  and  X£  is  also  gamma  (r,l)  truncated  to  the  right 

of  x  .  Thus  we  have 
o 

x 

EtXil  =  H  =  J  °  X  d^dx 


1 

F(xq) r (r) 


r  -x  , 
x  e  dx 


1 

(r) 


{r  I 


i=0 


r  1 

(r-i)  i 


r-i 


]} 


VAR [ X^ ]  =  a 2 


F  (x 
F(xo 


U1 


2 


X 


V(r)  o' 


r+l  -x  ,  2 

x  e  dx  -  u1 


1 

(r) 


{ (r+l) 


e 


-x 

o 


r+l 


t  l 

i=0 


(r+l) I 
(r+l-i) ! 


r+l-i 

o 


]} 


2 

“  yl 
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00 


e[X2]  = 


var[x2]  = 


Let 


A  = 


B  = 


F  (x)  -  F(x  ) 

1 2  *  !  X  d  1  -  F  (x  ) 

xo  ° 


— FT — T  /  xr  e  K  dx 

*2r(r)  x 

o 


1  r  Xo  f  rl  r-i, 

(e  l  tTTTTT  xn  1 


7T2r  i=0  (r_i^  1  ° 


2  F (x)  -  F(xq) 

rn^r- 


/  x  d  i  _ 


1  (  r+1  -x  j 

vtttJ  x  e  dx 

o 


1  -xo  — (£±111  x  r+l-i  2 

"2r(r)  6  iiO  (r+1_i)!  °  1,2 


r 

(r-i)  ! 


i 


x  r+l  i  i i \  i  ,i 

o  r  (r+l) !  r+l-i 
(r+l-i)  !  o 
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and  substitute,  then  we  have 


E[XX]  = 


r!  -  A 
(r )  ' 


VAR[XX] 


=  o 


2  =  (r+1) !  -  B  _  (rl  -  A) 


(r) 


2  2 

r(rr 


E[X2]  = 


TT2r(r)  ' 


VAR[X2] 


=  0 . 


B 


ir,r(r)  2  ,  .  2  * 
2  tt2  r(r) 


If  we  use  these  formulae  in  Theorem  2  in  Section  III  then 
we  have 


and 


U2r!  -  A)2 
2  (r  (r)  )  2r 


e  = 

ai  -  *1 

f 

ir2 

Thus 

(a^  ~  ^1^  ^u2r 1  "  A) 2 

(VI-A-1) 

P  = 

r  1T1  7t 2 2  (T  (r)  )  2 

For  given  p , 

x  and 
o 

r. 

ai  = 

2 

p  r!  (r-1)  !  n1  tt 2 

2  +  ^1 

(T72r!  -  A) 

(VI-A-2a) 
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and 


17 2  "  Ul(1  "  al5 


(VI-A-2b) 


The  value  of  (VI-A-2a)  and  (VI-A-2b)  are  determined  with 
given  p,  xq  and  r.  Thus  for  given  marginal  distribution  and 
correlation  coefficient  p  we  can  find  xo  ranges  as  before. 
For  simplicity,  we  will  consider  them  with  marginal  distri¬ 
bution  gamma  (2,1).  From  the  formula  of  (IV-0-2b) , 


*i  -  F<V 


-  X 


o 


e 


-x 

o 


IT 


2 


-x 


+  x 


e 


-x 

o 


and 


-x 


A  = 


o,  2  , 

(xo  + 


2x  + 
o 


2) 


Then  (VI-A-2a)  and  (VI-A-2b)  become 


a 


1 


a 


2 


2  p 


o 


+  ir . 


+ 


2  p 


(VI -A- 3a) 


(VI -A- 3b) 


These  and  are  probabilities,  so  they  have  to  be  greater 
than  or  equal  to  zero  and  less  than  or  equal  to  one; 
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2 


0 


< 


0 


-  2 


2  P  TT1  TT2 


+  IT. 


<  1 


*2  + 


2  p  ffl  ^2 

.  -2x 

4  o 

x  e 
o 


<  1 


From  these  two  inequality  equations  we  can  find  the  xq 
ranges  for  given  correlation  coefficient  p.  If  p  >_  0,  then 
and  a ^  are  always  positive.  Thus  we  only  need  to  find  the 
xQ  ranges  which  makes  (VI-A-3a)  and  (VI-A-3b)  less  than  one. 
From  equation  (VI-A-3a) ,  the  case,  and  from  equation 
(VI -A- 3b)  ,  the  a. ^  case,  we  have  exactly  the  same  inequality 
equation; 


x 

2pe  °(1 +  xQ) 


+  2p  (1 + 


xo} 


(VI-A-4) 


If  p  <_  0,  then  and  a are  always  less  than  one.  Thus  we 
only  need  to  find  the  xq  ranges  which  makes  (VI-A-3a)  and 
(VI-A-3b)  greater  than  zero.  From  the  equation,  we  have 


/27*(  1  +  xQ)  <_  xq2 


( VI-A-5a) 


where 


P*  =  |  P | 


and  from  the  a ^  equation,  we  have 
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2 


(VI-A-5b) 


_  x 

/2p  *  e  ° 


<  x 
—  o 


+  /2fT*(l  +  X  ) 
o 


where 

P*  =  |  P  | 

Note  that  if  we  let  for  the  left  side  equation  of  inequality 
sign  and  for  the  right  side  equation  of  inequality  equa¬ 
tions  (VI-A-4) ,  (VI-A-5a)  and  (VI-A-5b) ,  then  we  know  that  all 

y^  and  y^  a re  monotone  increasing  functions.  Thus  as  in  the 
exponential  case  we  can  find  the  range  of  xq  by  the  Newton 
Raphson  method.  For  the  positive  cirrelation  case,  equation 
(VI-A-4)  will  give  the  x^  and  and  for  the  negative  corre¬ 
lation  case,  equation  (VI-A-5a)  gives  x^  and  equation 
(VI-A-5b)  gives  x^.  By  using  subroutine  GBOUND  which  is 
listed  in  the  appendix,  we  can  find  x£  and  xu.  The  tables 
(Vl-a)  and  (Vl-b)  show  the  results  from  subroutine  GBOUND. 
Unfortunately  we  have  limitations  on  the  possible  range  of 
correlations;  we  can  generate  bivariate  random  vectors  with 
correlations  in  approximately  the  range  -.5  £  p  ^  .6  for  the 
Gamma  (2,1)  marginal  distribution.  In  the  same  way  as  the 
above,  the  values  x^  and  x^  can  be  computed  with  increasing 
complexity  for  integer  values  of  r  greater  than  2.  In 
principle  it  is  also  possible,  using  numerical  integration, 
to  compute  allowable  values  for  non-integer  values  of  r.  The 
computation  is  not  as  complicated  as  those  required  for 
Schmeiser's  bivariate  gamma.  Moreover  the  mixture-truncation 
method  does  not  require  that  the  user  be  able  to  compute 
the  inverse  gamma  distribution  function. 
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Table  Vl-a:  x  ranges  with  Gamma  (2,1) 
marginal  distribution 

and  given  p 


p 

XL 

X 

u 

P 

*L 

X 

u 

0.1 

0.41 

7.38 

-0.1 

0.93 

3.41 

0.2 

0.65 

6.23 

-0.2 

1.18 

2.75 

0.3 

0.89 

5.42 

-0.3 

1.35 

2.34 

0.4 

1.15 

4.74 

-0.4 

1.50 

2.03 

0.5 

1.47 

4.09 

-0.5 

1.62 

1.79 

0.6 

1.95 

3.32 

Table  Vl-b:  x  ranges  with  Gamma  (3,1) 
marginal  distribution 

and  given  p 


P 

*L 

X 

u 

P 

XL 

X 

u 

0.1 

0.82 

8.99 

-0.1 

1.64 

4.73 

0.2 

1.17 

7.72 

-0.2 

1.37 

3.98 

0.3 

1.49 

6.82 

-0.3 

2.21 

3.5 

0.4 

1.83 

6.05 

-0.4 

2.4 

3.15 

0.5 

2.24 

5.31 

-0.5 

2.56 

2.89 

0.6 

2.84 

4.42 

135 


B.  GENERATING  PROCEDURE 


The  generating  procedure  of  bivariate  gamma  random  vec¬ 
tors  is  almost  the  same  as  the  uniform  and  exponential  cases. 
We  also  developed  here  three  methods  which  are  the  FXO 
method,  the  UXO  method  and  the  TXO  method.  As  in  the  uniform 

and  exponential  cases,  these  three  methods  are  the  same 

except  in  the  choice  of  x  from  the  x  range  [x.,x  ].  The 
FXO  method  chooses  xq  as  a  fixed  point  from  the  xq  range  and 
uses  it  during  the  entire  routine  while  the  UXO  and  TXO 
methods  choose  another  xq  in  every  routine.  From  experience, 
the  midpoint  of  the  xq  range  gives  the  best  result  of  the 

FXO  method.  In  the  UXO  method  and  the  TXO  method,  we  assume 

that  xq  has  the  uniform  distribution  and  triangular  distri¬ 
bution,  respectively,  over  [x„,xu].  This  assumption  removes 
some  discontinuity  which  occurs  at  the  truncation  point  in 
the  FXO  method. 

Gamma  Mixture-Truncation  Method 
(for  integer  shape  parameters) 

1.  (Initialization) 

i)  For  given  allowable  correlation  coefficient  p. 


find  x.  and  x 

£  ■ 


u 


2.  Define  truncation  point  x 


o 


*  FXO  method 


i)  x 


o 
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*  UXO  method 


i)  Generate  a  uniform  [0,1]  random  variable 
ii)  Xo  =  Xjl  +  (xu-xe)  * 

*  TXO  method 

i)  Generate  two  uniform  [0,1]  random  variables, 

V  V2 


ii)  x  =  x  +  x.  +  x» 
o  £  1  2 


where 


x  =  (x  +  x  )  /2 . 0  , 
m  a  u  ' 


X1  =  ■ 


x-  =  (x  -  x  ) *V~ 

2  u  m  2  * 


3.  Compute  parameter  values 


-x 


r-1  e  x 


o  j 


-  F<V  =  1  -  l 

j=o 


I  r 


=  1  - 


p  r !  (r-1)  !  2 


(^r!  -  A) 


2  +  ^1  ' 


77 1 

1  -  ^(1  -  «i>  ■ 


where 
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A 


e 


r 


r-i 


-x 


r 

l 

i=0 


(r-i) 


x 


4.  Choose  type  for  Y 

i)  Generate  a  uniform  [0,1]  random  variable  U 

ii)  If  U  <  go  to  9 

5.  Y  is  an 

i)  Generate  a  gamma (r,l)  random  variable  G^ 

ii)  If  G^  >  xq/  set  Y  G^  and  go  to  6 

iii)  Otherwise  return  to  5.i) 

6.  Choose  type  for  Z 

i)  Set  U  ^  (  (U  -  771)/(1  -  t^)  ) 

ii)  If  U  <_  1  -  »2#  go  to  8 

7.  Z  is  an 

i)  Generate  a  gamma  (r,l)  random  variable  G2 

ii)  If  G2  >  xq,  set  Z  +■  G2  and  go  to  11 

iii)  Otherwise  return  to  7.i) 

8.  Z  is  an 

i)  Generate  a  gamma  (r,l)  random  variable  G2 
ii)  If  G.  <  x  ,  set  Z  +■  G„  and  go  to  11 

w  O  4 

iii)  Otherwise  return  to  8.i) 
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9 .  Y  is  an 

i)  Generate  a  gamma  (r,l)  random  variable  G^ 

ii)  If  G^  <_  xq/  set  Y  G^  and  go  to  10 

iii)  Otherwise  return  to  9.i) 

10.  Choose  type  for  Z 

i)  Set  U  U/i^ 

ii)  If  U  <_  a^,  go  to  8 

iii)  Otherwise  go  to  7 

11.  Deliver  (Y,Z)  and  go  to  4  for  the  FXO  method,  or  go 
to  2  for  the  UXO  and  TXO  methods  until  a  sufficient 
number  of  random  vectors  are  obtained. 

For  the  non-integer  shape  parameter  case,  step  1  and  step 
3  need  more  complicated  computation  procedures  but  here  we 
are  only  concerned  for  integer  shape  parameter  cases.  In 
step  1  we  used  the  subroutine  GBOUND  which  is  listed  in  the 
appendix.  The  scatter  plots  and  bivariate  histograms  of  the 
resulting  random  vectors  from  this  algorithm  are  shown  in 
the  next  section. 

C.  SIMULATION  RESULTS 

The  scatter  plot  and  bivariate  histograms  for  random 
vectors  of  the  mixture-truncation  bivariate  gamma  variables 
with  marginal  gamma  (2,1)  distribution  and  given  correlation 
(0.3  and  -0.3)  are  given  here.  As  in  the  uniform  and 
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exponential  cases,  the  UXO,  TXO  and  FXO  methods  are  used 
and  are  shown  in  Figures  (Vl-a)  ,  (Vl-b)  and  (VI-c)  ,  respectively. 

The  FXO  method  has  some  discontinuity  at  the  truncation 
point,  but  the  UXO  and  TXO  methods  appear  to  give  a  relatively 
smooth  continuous  distribution.  The  computed  correlation  in 
subroutine  BIVHST  (Bivariate  histogram)  is  a  little  differ¬ 
ent  from  the  given  correlation.  But  we  can  assume  this 
difference  is  a  result  of  sampling  error. 


140 


Figure  Vl-al. 


Scatter  plot  for  sample  of  size  2000 

from  mixture-truncation  bivariate  Gamma 

distribution  with  p  =  0.3.  Here  x  has 

uniform  distribution  over  (xn/x  ).  ° 

a  u 
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Figure  VI-a2. 


Bivariate  histogram  for  sample  of  size  2000 
from  mixture-truncation  bivariate  Gamma 
distribution  with  p  =  0.3. 
uniform  distribution  over 


Here  x  has 


<Vxu>' 
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Figure  VI-a3.  Scatter  plot  for  sample  of  size  2000  from 

mixture- truncation  bivariate  Gamma 

distribution  with  p  =  -0.3.  Here  x  has 

uniform  distribution  over  (x  ,x  ).  ° 

a  u 
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Figure  VI-a4 .  Bivariate  histogram  for  sample  of  size 

2000  from  mixture-truncation  bivariate 

Gamma  distribution  with  p  =  -0.3.  Here 

x  has  uniform  distribution  over  (x  ,x  ). 
o  i'  u 
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Figure  Vl-bl.  Scatter  plot  for  sample  of  size  2000 

from  mixture  truncation  bivariate  Gamma 
distribution  with  p  =  0.3.  Here  x  has 
triangular  distribution  over  (x^ ,x°)  . 
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Figure  VI-b2 . 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
Gamma  distribution  with  p  =  0.3.  Here 
x  has  triangular  distribution  over 

(xt'xuK 
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Figure  VI -b3.  Scatter  plot  for  sample  of  size  2000  from 

mixture-truncation  bivariate  Gamma 
distribution  with  p  =  -0.3.  Here  x  has 
triangular  distribution  over  (x?/xuJ. 
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Figure  VI-b4. 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
Gamma  distribution  with  p  =  -0.3.  Here 
xq  has  triangular  distribution  over 
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Figure  VI-cl.  Scatter  plot  for  sample  of  size  2000  from 

mixture-truncation  bivariate  Gamma 

distribution  with  p  =  0.3.  Here  xQ  is 

fixed  at  the  midpoint  of  x  and  x  . 

£  u 


149 


0002  •  ins  IldHVS  AtlVtaHnS  I1<!4¥j  HVlbVMO 


M  U> 


U  u. 
CO 

— c 

O.W 


4K4 

>0.0. 

OCA 


r  ■* 
at  *- 
Ovo 


MOO 
aw  >*• 
— U'uM 
MJO 


? 


►  • 
wX» 


uoo 

■4  -0-4 


OO  OOO  OO  Oai 
OO  OOO  Oo  OO 
» 

U  U  Ui  U  U*U  u  u  oju- 
_,V  O  ***«* 

a  n  O  Oc'fl  oc  oa 
y  or* 

«-  i*io  ^rrvjfv,  <mi—  «**'# 

vi-c-o  w  -o>r  oa-  av< 
C“  -C  «*  o  a>  rr  C'  «V*l 
>••  •  •  •  »  *  •• 
•*<«  —  (M  «h# 


2 


u  ►■ 


M  «_>*••  OLlV*  uu  cy 

►  CO  oou  CJO  oo 


u  u  UI 

-**»  o 

o  a:  O 
X  -  c 


X  •  • 

r>*< 


uju  u  u<u 
r-t/\f*»  r~vn 

IVN-I  — o 
K'O'f  O  J- 

O  —INI  (fO 

k><00  «N*f* 

CM-O 


U  tt» 
-O 


OO 

CO 

-re 

INfO 


4C11I  20 

—  O  5K 
aC2  u'c£ 


a  i 

20 

XX 


Uiu.  <i»-  «J  *3  <»- 
XX  ><^>0t  i/»X  XX 


I  I  I 

: ; s 

!! 


s: 


ri 


\  t 


IMI 
lilt 
l  I  I  l 
till 


III 

III 

I  I  a  ■ 


I  I  I  I  I  I  I  I 


I  I  I  I  l  I  I  I 


i  I  l 
I  I  I 


I  l  I  I  I  i  I  l  I 


I  I  I  I  I  l  •  l  l 


till 
I  I  l  I  «  I 

I  l  I  I  I  I 


l  I  I  I 
l  I  I  I 
I  l  I  I 
I  I  l  I 


i; 


Figure  VI-c2. 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
Gamma  distribution  with  p  =  0.3.  Here 

x  is  fixed  at  the  midpoint  of  x^  and 

O  36 


X 


u 


150 


* 

X 


* 


Figure  VI-c3. 


Scatter  plot  for  sample  of  size  2000 
mixture  truncation  bivariate  Gamma 
distribution  with  p  =  -0.3.  Here 
fixed  at  the  midpoint  of  x,  and 
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Figure  VI-c4. 


Bivariate  histogram  for  sample  of  size 
2000  from  mixture-truncation  bivariate 
Gamma  distribution  with  p  =  -0.3.  Here 
is  fixed  at  the  midpoint  of  x  and  x 
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NO.  CeSERVATI?N$ 


VII.  CONCLUSION 


The  mixture- truncation  method  is  a  general  method  which 
can  generate  bivariate  random  vectors  having  any  theoretical 
marginal  distribution  and  allowable  correlation.  The 
generating  procedure  is  very  simple  and  doesn't  need  much 
computation  for  defining  parameter  values.  In  this  respect, 
the  mixture-truncation  method  is  a  very  attractive  method  for 
generating  bivariate  random  vectors.  A  price  is  paid  for 
this  simplicity  and  generality  in  that  the  Frechet  bounds  of 
correlation  for  the  bivariate  distributions  specified  by  the 
marginal  distribution  given  by  Moran  (1967)  are  not  always 
attained.  Also  there  is  some  discontinuity  in  the  bivariate 
distribution.  However  this  discontinuity  can  be  decreased 
by  giving  some  distribution  to  the  truncation  point  over 
its  range  for  given  p .  Thus  the  mixture-truncation  method 
is  very  attractive  for  simulation  studies  involving  only 
partly  specified  dependency  structures.  The  mixture- 
truncation  method  may  be  extended  to  generate  bivariate 
random  vectors  having  negative  values.  Another  extension 
may  be  made  to  use  grade  correlation  or  rank  correlation 
which  are  invariant  under  transformation  instead  of  using 
the  product  moment  correlation. 
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APPENDIX:  PROGRAM  LISTINGS 


THE  PROGRAM  FOR  GENERATING  9 1  VAR  I  AT  E  RAOQM  V=C T CRS < Y, Z ) 
HAVING  IDENTICAL  MARGINAL  DISTRIBUTION  BY  M I XT U R E-TRUN C A 
C-TICN  METHOD. 

THE  CORRELATION  LIMItc  ARE  AS  FOLLOWS  FOR  GIVEN  MARGINAL 
DISTRIBUTION. 

1.  UNIFORM:  -0.75  <  CR  <  0.75 

2.  EXPONENTIAL  :  -0.48  <  CP  <  0.64 
i .GAMMA  (2,1) :  -.55  <  CR  <  0.64 

GENERATING  PROCEDURE 

CALL  BUNF  ( CR  ♦  NV  ,Y  ,  Z  ,  ! S  ) 

CALL  B5XP(CR, NV,Y, Z , IS ) 

CALL  BGAM2(  CR  »NV  , Y, Z»I S) 

WHERE 

CR:  GIVEN  CORRELATION. 

NV:  NUMBER  OF  RANDOM  VECTORS  BE  GENERATING. 

Y:  THE  FIRST  ELEMENT  OF  VECTOR  WITH  DIMENSION  NV. 

Z:  THE  SECOND  ELEMENT  CF  VECTOR  WITH  DIMENSION  NV . 

IS:  INITIAL  SEED  PQR  UNIVARIATE  GENERATOR. 


SUBROUTINE  BUN F ( CR , NV , Y , Z , I U ) 

C 

C 

DIMENSION  Y  (  N  V  ) » Z  (  N  V  ) 
r 

IF(CR  .GT.C.75  .OR.  CR.LT.-Q.75)  GC  TC  70 

Ic  (CR.LT.0.0 )  GO  T0  5 

XL  =  C • 5-< (  (9.0-12.0*CR)**0.5)/6.0) 

XU=  1 .  O-XL 
GC  TO  10 

5  XL  =  (  ( -1 .  0*C  R )  /3  .0  )  **0 .5 
XU= 1 • O-XL 
10  CO  100  1=1, NV 

call  random < i l,u, i  ) 

X0=  XL  +  ( XU  — XL  )  *U 

CALL  UP  JO  B(  CR  ,X0,  A  PI,  AP2,  PI  E1,PI  E2,  BETA  ) 

C - CHOOSE  TYPE  FOR  Y 

CALL  RANDOM  (  IU,U,T  ) 

IF(L  .LE.PIS1  )  GO  TO  50 
L= ( U-PI E 1 ) / (  1.0-PIE1) 

AP=  1 .  O-AP  2 
C  Y  PRO N  X2 

20  CALL  RANDOM  (  !  L  ,W  ) 

Y ( I )=XO+W*( l.C-XO) 

C  CHCJSE  TYPE  FOR  Z 

IF(U.LE.AP)  GC  TO  40 
C...Z  FROM  X2 

30  CALL  RANDCM( I L, W, 1 ) 

2(1  )=X0+W*(1.C-X0) 

GO  TC  100 
C...Z  FROM  XI 

40  CALL  RAN  DOM (  IL,W, 1) 

Z(I )=XO*W 
GO  TC  100 
C... .Y  FRCM  XI 
50  U  =  U  /  PI E 1 

60  CALL  RANOOM ( IU,W ,1  ) 

Y(I  )=XO*W 

IF(U  .LE.AP1  )  GO  TO  40 
GO  70  30 
100  CONTINUE 
RETURN 

70  WR  ITE(6 ,8  C)  CR 

80  F0RMAT(5X, 'CHECK  CORRELATION  LIMIT, CO  NOT  ALLOW  GIVEN, 
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CCTPREIAT ICN  IN  THIS  °RCGRAM  FOR  CR=»,P7.3) 

STOP 

ENC 

C 

C 

SUBROUTINE  UP JOB ( C R , XO , AP 1 , AP 2 , P I E I , P I  52 , BET A  ) 
C - UNIFORM 

C - FIND  INITIAL  VALUE 

£°1=<CR  +  3.0:*XC**2)/(3.05*X0) 

AP2=1.0-< 1.0- API)* (XO/ (1.0- XO ) ) 

P I  E 1= ( 1.0-AP2 )/( 1.0-AP1+1.0-AP2) 
PIE2=(1.0-AP1)/<1.0-AP1+1  .0-AP2  ) 
eETA=APl-( 1 .0-AP2  ) 

RETURN 
ENC 


SUBROUTINE  8EXP(CR, NV,Y,Z,IS) 

C 

C 

DIMENSION  Y ( NV ) tZ(NV) 

C 

IX= I S*0 • 8 

IU= IS*0. 3+123456 

IF <CR .GT.0.64  .OR.  CR.LT.-0.48)  GC  TO  70 
CALL  BOUNC(CR  ,XS1 »XS2,XL,XU) 

CO  IOC  1=1, NV 
CALL  RANDOM ( IU ,U , 1  ) 

XO  =  XL  +  ( XU-XL) *U 

CALL  EPJOe<  CR,X0,AD1,AP2,PIE1,PIE2,E6TA) 

C - CFCOSE  TYPE  FOR  Y 

CALL  RANDOM ( IU,U, 1  ) 

IF(L.LE.PTEl)  GO  TO  50 
L= ( L-°  I E 1 )/( 1  .0-PIE1) 

AP= 1 . 0-AP2 

C - Y  FROM  X2 

2 C  CALL  EXPONdX  ,X,1  ) 

V< I )=X+XO 

C - CFCOSE  TYPE  FOR  Z 

IF(L.LE.AP)  GC  TO  4C 

C - Z  FROM  X2 

30  CALL  EXPONdX, X,l) 

Z(I )=X+XO 
GC  TC  100 

C - Z  FROM  XI 

40  CALL  RANDOM ( I L ,U, 1) 

2(1 )=-1.0*AL0G(1.0-U*°IEl  ) 

GC  TC  100 

C - Y  FROM  XI 

50  U=U/PI 51 

60  CALL  RANOCM (  IL,U,1  ) 

Y(I ) 1 • 0  *A  LOG ( 1. 0-U*P I E 1 ) 

IF(L.LE.AP1 )  GC  TO  40 
GO  TC  30 
100  CONTINUE 
RETURN 

70  ViR  I  T  E  (6 , 8  C )  CR 

80  F0RMAT(5X, ‘CHECK  CORRELATION  LIMIT, CC  NOT  ALLCW  GIVEN, 
CCOPPELATICN  IN  THIS  PROGRAM  FCP  CP=‘,P7.3) 

STCP 

END 

C 

C 

SUBROUTINE  EPJQB(CR,XO, API , AP 2,P I  El ,pI E2, BETA  ) 

C  EXPONENTIAL 

PIE1=1.0-(1.0/EXP(XC)) 

P I  E  2  =1 . 0-  PI  El 
AP1=PIE1*< 1 .0+(CR/XC**2  ) ) 

AP2=1.0-( (P IE  1/PTE 2  )*( 1 .0-AP1  )  ) 
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BcTfl*APl-( 1.0-AP2) 

RETURN 

ENC 


SUBROUTINE  BOUND <  RO ,X$ltXS2,XL»XU) 
IF(PC.LT.O.O)  GO  TO  10 
AA  = ( 1.0-RC- ( ( RO**  2 )* (5. 0/12.0))  )**0.5 
XS1=< ( 1 .O-Q.S^ROJ-AA )/(R0/3.0) 

XS2  =  ( (1  . 0-0.5 *R0) +  A A) /(R0/3.0 ) 

CALL  RPS0N( XS 1tXX,RC) 

X  l=  X  X 

CALL  RPSON(XS2,XXtRC) 

XU=XX 

RETURN 

10  AR  =  ABS( RO  ) 

XL= AR**0 . 5 
X$l=C.O 

XS2= ( ( AR**0 • 5 )— AR ) / ( AR*0.  5 ) 

CALL  RPS0N(XS2tXX,RC) 

XU  =  XX 

RETURN 

END 


SUBROUTINE  R P SON ( XS , XX , PO ) 

EE=C  .0001 
V 1  =  X  S 
N  =  1 

10  IF( RC.LT.O.O)  GO  TO  20 
CALL  FUNP(Y1»FV,DV,PC) 

GO  TO  25 

20  CALL  FUNN<Y1,FV,9V»R0) 

25  >2=  V  1-(FV/DV) 

IF |  (ABS( Y1-Y2  )) .LE . EE)  GO  TO  30 
Y1  =  Y2 
GO  TC  10 
30  XX=Y2 
RETURN 
END 
C 

c 

SUBROUTINE  FUNP(Y1»FV.DV»RQ) 

FV= (Y 1**2 )-(RC*( EXP ( Y1 ) -1 -0 ) ) 
CV=(2.0*Y1)-(R0*EXP(Y1  )  ) 

RETURN 

END 

C 

C 

SUBROUTINE  FUNN( Y1 , FV, DV, RO J 
RR  =  A  E S ( RO  ) 

FV  =  (Y1**2 )-PR*<  (EXP  ( Y 1  )  - 1 . 0  )  **2  ) 

CV=  2 • 0*Yl-( 2. C*RR*E  XP ( 2  .0*Y1)  )  +  ( 2  .0*RR *EXP ( Y 1  )  ) 
RETURN 
END 


SUBROUTINE  BGAM2 < C R  ,NV tY , Z , I S ) 

C 

C 

DIMENSION  Y  ( N V ) ,Z(NV) 

C 

IX=IS*0.8 
IU=IS*0. 3  +  123456 

IF(CR.GT.0.64  .OR.  CR.LT.-0.55)  GO  TO  70 
CALL  GF?0UND(CR»XS1,XS2tXL,XU) 

CO  100  I =1 » NV 
CALL  RANDOM ( I U  »U  ♦  1 ) 

XO=XL+( XU-XL) *U 
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CALL  GPJ08 (CP  ,XQ,AP1,AP2,PI51,PI E2*EETA) 

C - CFCOSE  TYPE  FOR  Y 

CALL  RANDCM< !L,U,1  ) 

IP(l.LE.PIEl)  GO  to  50 
U=(L-PIE1J/(1.0-PIE1) 

AP  =  1  .O-AP? 

C - Y  FRCNX2 

20  CALL  GAMA(2.C,IX,X,1) 

IF(X.LE.XC)  GC  TO  2C 
Y( I  )  =X 

C - CFCOSE  TYPE  FOR  Z 

IF(U.LE.AP)  GC  TO  40 

C - Z  FROM  X2 

30  CALL  GAMA(2.0,IX,X,1) 

IF(X.LE.XC)  GC  TO  3C 
Z ( I ) -X 
GO  TQ  100 

C - Z  FRCN  XI 

40  CALL  GAMA(2.0,IX,X,1) 

IF(X.GT.XC)  GC  TO  4C 
Z( I  )  =  x 
GO  TO  100 

C - Y  FRO  N  XI 

50  L=U/PI El 

60  CALL  GAMA(2.0,IX,X,1) 

IF(X.GT.XO)  GC  TO  6C 
Y  (  I  )  =  X 

IF<L.LE.AP1 )  GO  TO  40 
C-0  TC  30 
100  CONTINUE 
RETLRN 

70  VnRI7E(6,80)CR 

80  FORMAT (5X  » 1  CHECK  CORRELATION  LIMIT, CC  N0T  ALLCV,  GIVEN, 
CCORRELATICN  IN  THIS  PROGRAM  FOR  CR=',F7.3) 

STOP 

END 

C 

C 

SUBROUTINE  GP JOB ( CR , XO , AP 1 , AP 2 , P I E 1 , P 1 52 , BE ta ) 

C  GAPMA(2»X) 

EX=1.0/EXP(X0) 

F IE  1  =  1 .0-EX-XC*EX 
PIE2=l.0-PIEl 

AD=(X0**4)*(1 .0/EXP(2.0*XQ)  ) 

API =PI E1+ ( (  2.C*CR*PIE1*PIE2**2)/AC) 

AP2=1.0-( (PIE  1/PI F2  )*( 1  .0-AP1 ) ) 

BETA=APl-( 1 .0-AP2 ) 

FETLRN 

END 

C 

C 

SUBROUTINE  GB CUND  (  RC  ,X  S  1  ,XS 2 ,  XL ,  XL'  ) 

IF(RO.LT.O.O)  GO  TO  10 

A1 = ( (4. 0* PO ) / 3 .0 ) +SQRT (  (<4.0*R0**2  ) / 9 .0 )  +  ( 4 . 0 *RC ) ) 
A2=2.0*(1.0-(RO/3.0)) 

XS1=A1/A2 
XS1=XS1+1  .0 
>S2=XS1+1 C.O 
CALL  GRP$0N(X$1,XX»PQ) 

XL  =  X  X 

CALL  GRPSON( XS2,  XX, RO) 

XU  =  XX 
RETURN 

10  AR=2.0*ABS(RC) 

XL=($GRT(  AR)+$QRT( AP+SQRT( AR ) *4.0  )  )/2.0 
XS1=C.0 

A1=SQRT( (SQRT(AR)/6.0>+(R0/9.0) )- (SCRT( AR ) /6. C) 

A2=SGRT(AR)/12.0 

X$2=A1/A2 

CALL  GRPS0N(XS2,XX,P0) 

XU  =  X  X 
RETURN 
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END 


C 

C 

SU8 PCUTINE  GRP  SO  NOS, XX,  RO) 

EE=C .001 
V  1  =  X  S 

10  IF (RC.LT.0.0)  GO  TO  20 
CALI  GFUNP(Y1  ,FV,DV,RQ) 

C-C  TO  25 

20  CALL  GFUNN(Y1  ,FV,DV,PC) 

25  >2=  U-(FV/DV) 

IF< (ABSIY1-Y2  ))  .LE.EE)  GO  TO  30 
Y1  =  Y2 
GO  TO  10 
30  XX=Y2 
RETURN 
END 
C 
C 

SUBROUTINE  GFLNP( Y1  ,FV,DV,RO) 

Fl  = ( Yl**4)  +  (2 .0*RO* ( 1.0+Y1 )**2.0 ) 
F2=2.0*R0*EXP(Y1)*( 1  .O  +  Yl ) 

F  V  =  F 1-F2 

CF1=(4.0*Y1**3)+(4.C*R0*( 1.0+Y1) ) 
CF2=2.0*RC*EXP(Y1)*(2.0+Y1) 

CV  =  C  F 1-DF  2 
RETURN 
END 
C 
C 

SUEFCUTINE  GFUNN (Y1  ,FV, CV,RO) 
TR=SQRT(2.0*A3S(R0)  ) 
F1=<Y1**2.0)+(TR*( 1.0+Y1) ) 
F2=TR*EXP(Y1  ) 

FV  =  F  1-F  2 

CF  1  =  ( 2 . 0*Y1 i  +  TR 

CF2=TR*EXP( Y1  ) 

CV=  C  F 1-DF  2 
RETURN 
END 


SUBROUTINE  BIVHST  (X,  Y,  N,  WORK) 

C 

IMPLICIT  REAL*8  (D) 

INTEGER*2  CELL,  MAX,  SCALE,  ONE,  CNE5 
INTEGER  RUNS,  WN,  SN,  WEIGHT,  U,  FLAG 

REAL  D,  DELTA,  DELTAX 

C 

DIMENSION  X  ( N  )  ,  Y ( N  )  ,  WCRK(N) 

DIMENSION  XM0M6),  YM0M(6),  OWCRK(E),  CELL(32,32) 
DIMENSION  KEY (16) , XLABEL! 8) 

DATA  NOUT/6/,  ONE/1/,  0NE5/15/ 

C 

IF(N  . GT .  15)  GO  to  20 
WRITE! NOUT , 10  )  N 

10  FORMAT  (  'l^  +  BIVHST***  TOO  FEW  SAMPLE  POINTS.  N  =',IS) 
RETURN 
C 

20  AN  =  N 

AN2  =  AN  *  AN 
NM  =  N  -  1 

C - CRDER  X-SAMPLE  AND  FIND  RANGES 

CALL  SORTON  0,1, N,Y) 

XR ANGE  =  X<  N)  -  X(1  ) 

YM  A  >  =  Y (  1) 

YM  IN  =  Y (  1 ) 

CO  30  1=2, N 

I F ( Y (  I  )  .GT.  YMAX )  YMAX  =  Y(I) 
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I F ( Y  ( I  )  .IT.  YMIN)  YMIN  =  Y(I) 

30  CONTINUE 

YRANGE  «  YM AX  -  YMIN 

IP ( XRANGE  *  YRANGE  .GT.  0.)  GO  T0  50 
WRITE (NQUT , 40 ) 

40  FORMAT ( ,1***B!VHST***  X  OR  Y  SAMPLE  HAS  ZERO  RANGE' ) 
RETURN 
C 

C - ZERC  OUT  CATA  ARRAY 

5C  DO  7 C  1=1  ,32 
CO  6C  J  =  1  ,32 
CEL  L ( I  ,  J )  =  0 
60  CONTINUE 
70  CONTINUE 

C - FINC  SCALE  FACTORS 

AX  =  31.999  /  XRANGE 
EX  =  1.0001  -  AX  *  X(1  ) 

CELTAX  =  XRANGE  /  32. 

AY  =  31.999  /  YRANGE 
PY  =  1.0001  -  AY  *  YMIN 
CELIA  =  YRANGE  /  32. 

MAX  =  0 

C - ACCUMULATE  CELL  COUNTS 

DO  100  1=1, N 
IX  =  AX  *  X( I  )  +  BX 

IY  =  AY  *  Y  C I  I  +  BY 

IY  =  33  -  IY 

CELLCIX,  IY)  =  CELL (  IX,  IY)  +  ONE 
IF(CELL(IX,  IY)  .GT.  MAX)  MAX  =  CELL( IX,  IY) 

100  CONTINUE 

C - SCALE  FACTOR  FOR  COUNTS 

SCALE  =  (MAX  -  ONE)  /  0NE5  +  ONE 

C - FINC  SAMPLE  MOMENT S 

CWORK(l)  =  O.C 

CALL  ACCUM(X,N,DWORK) 

CALL  MOMENT (DWORK,  XMOM) 

OWORK(l)  =  O.C 

CALL  ACCUM(Y,N,DWORK) 

CALL  MOMENT (DWORK,  YMOM ) 

C - FINC  COVARIANCE 

CSUM  =0.0 
DO  110  1=1, N 

CSUM  =  OSUM  +  ( X ( I  )  -  XMOM(l))  *  <  Y  C  I )  -  YMOM(l)) 

110  CONTINUE 

COVAR  =  DSUM  /  AN 
XDEV  =  SORT ( XMQM( 2)  ) 

YOEV  =  SORT ( YNOM<  2 )  ) 

RHO  =  COVAR  /  (XOEV  *  YDEV) 

C - KENCALL'S  TAU  COEFFICIENT 

T  =  C.O 

IF(N  .GT.  500)  GO  TC  119 
CO  118  1=1, NM 
00  114  J=  I  ,  N 
2  =  Y(J)  -  Y  C I ) 

IF ( Z  .GT.  0.)  T  =  T  +  2.0 
I F  ( Z  -LT.  0.)  T  =  T  -  2.0 
114  CONTINUE 

118  CONTINUE 

T  =  t  /  ( AN2  -  AN) 

C - SET  UP  WORK  AREA  AND  ORDER  Y  SAMPLE 

119  CONTINUE 
VI  =  1.0 

CO  120  1=1, N 
WORK  (I)  =  WI 
W I  =  WI  +  1.0 

120  CONTINUE 

CALL  SORT  ON ( Y ,  1,  N,  WORK) 

C - FINC  RANK-CRDER  STATISTICS 

VN  =  0 
SN  =  C 
RUNS  =  0 

I F (  Y(  1)  .LT.  XI )  )  RUNS  =  1 


159 


130 


140 


AN 


N)  GO  TO  200 


150 


SN  = 
IX  = 


160 


C  170 


180 


190 


FLAG  =  0 
C  =  C.O 
!  =  1 
IX  =  1 
IY  =  1 

Z  =  ABSIFLOAT (IX-TY  )  )  / 

IF ( C  .IT.  Z)  C  *  Z 
I F ( X ( IX  J  -  Y (  IY )  )  150,  160,  140 

•CURRENT  VALUE  IN  POCLEO  SAMPLE  IS  A  Y 
RUNS  =  RUNS  +  FLAG 
FLAG  =  0 
IY  =  IY  +  1 
1=1  +  1 
I F  <  I Y  .GT 
GG  TC  130 

•CURRENT  VALUE  IN  POCLEO  SAMPLE  IS  AN  X 
RUNS  =  RUNS  +  1  -  FLAG 
FLAG  =  1 
WN  =  WN  +  I 

SN  +  WEIGHT ( I , N ) 

IX  +  1 

IF ( I X  .GT.  N)  GO  TO  220 
1  =  1  +  1 
GG  TG  130 

•X  AND  Y  VALUES  ARE  TIED.  USE  MIDRANK  CORRECTICN 
BASE  =  X( IX) 

N  =  1 

J1  =  I  +  I  +  1 

J2  =  WE IGHT( I  ,N )  +  WEIGHT ( I +1  ,N ) 

K  =  I 
1=1+1 

•GET  TIED  X  VALUES 
IX  =  IX  +  1 

!F(  IX  .GT .  N)  GO  TO  180 
IF ( X ( I X )  .NE.  BASE)  GG  TO  180 
N  =  N  +  1 
1  =  1  +  1 
J1  =  J1  + 

J2  =  J2  + 

GG  TC  170 

■FINC  TIED  Y  VALUES 
IY  =  IY  +  1 

IF (  IY  .GT.  N) GO  TO  190 
I F  <  V  <  IY )  .NE.  BASE)  GO  TO  190 
1  =  1  +  1 
J1  =  J1  +  I 
J2  =  J2  +  WEIGHT! I  ,N) 

GO  TO  180 

•WEIGHTED  AVERAGE  OF  TIED  VALUES 
L  =  I  -  K  +  1 

Jl)  /  FLCAT(L)  +  0.5 


WEIGHT ( ! ,N) 


J2)  /  FLOAT (L )  +  0.5 


( L-M)  )  /  =LOAT(L) 


Z  =  FLO  AT (  M  * 

WN  =  WN  +  Z 
Z  =  FLOAT ( M  * 

SN  =  SN  +  Z 

- AC  FCC  CORRECTICN  FCR  RUNS 

RUNS  =  PUNS  +  FLOAT ( 2  *  M  * 

1  =  1+1 

IF ( I X  .GT.  N)  GO  TO  22  C 
IF (  IY  .LE.  N)  GO  T0  130 

- Y  CESERVATICNS  EXHAUSTED 

200  RUNS  =  RUNS  +  1 
N2  =  N  +  N 
DC  210  1=1, N2 
WN  =  WN  +  I 
SN  =  SN  +  WEIGHT ( I , N  ) 

210  CONTINUE 

- FINC  MANN-WHITNEY  STATISTIC 

22C  L  =  WN  —  (N  *  (N  +  1))  /  2 

- SPEARMAN'S  RANK  CORRELATION  COEFFICIENT 

CSUN  =  0. 

WI  =  1.0 
DO  230  1=1, N 


160 


oooo 


CSUM  =  DSUM  +  (WI  -  WORK ( I ) )  **  2 
WI  =  WI  +  1.0 
230  CONTINUE 

R  =  1 .  -  6.  *  DSUM  /  (AN  *  ( A  N2  -  1.)) 

C - FINC  SAMPLE  MEDIANS 

I F  (  ROD ( N  *  2 )  .EQ.  0)  GO  TO  240 
M  =  1  +  N/2 
XMEC  =  X (  M ) 

YMEC  =  Y  (  M  ) 

GO  to  250 
240  M  =  N/2 

XMEC  =0.5*  ( X ( M |  +  X<  M+l)  ) 

YMEC  =  0.5  *  (Y(M)  +  Y( M+l) ) 

C - FINC  NORMALIZED  STATISTICS 

250  RNCPM  =  (RUNS  -  N  -  1)  *  SORT ( ( AN  +  AN  -  1.)  /  ( AN2  - 
CAN)  ) 

FAC  =  SORT ( 12 .  /  ( AN2  *  (AN  +  AN  +  1.))) 

UNORM  =  (U  -  C.5  *  AN2)  *  FAC 
VNOPM  =  (WN  -  AN2  -  0.5  *  AN)  *  FAC 

SNCPM  =  ( SN  -  AN2  -  0.5  *  AN)  *  FAC 

C - REORDER  Y  SAMPLE 

CALL  SORTON  (WORK,  1,  N,  Y) 

WRITE  OUTPUT 

- HEADING 

WRI TE (NOUT ,  3C0 )  N 

300  FORMAT!  •  1  »  ,20X,' BIVARI  ATE  SAMPLE  SUMMARY ', 10X ,' SAMPLE 
C  SIZE  =  ',  19//  ) 

WRITE!  NOUT  ,310 

310  FORMAT ( 13 X, •+--• ,7( ,7('-'l ), •* - +  •  ) 

C 

C - ECDY  OF  PLOT 

CALL  OUTPUT ( C  ELL ( 1  ,  1),  SCALE,  1,  C.) 

CALL  REPEAT 

CALL  QUTP  L'T  (  C  ELL  (  1  ,  2),  SCALE,  2,  Y  M  AX  -  DELTA) 

WR ITE ( NOU  T , 32  C ) 

320  FORMAT ('+', 93X ,' MEASURES  OF  ASSOCIATION') 

CALL  REPEAT 

CALL  OUTPUT ( C  ELL(  1 ,  3),  SCALE,  3,  C.) 

WRITE! NOUT ,330)  COVAR 

330  FORMAT(  •  +  • ,83X, 'COVARIANCE' ,2 5X, 1PE  14.6) 

CALL  REPEAT 

V* R I  T E  <  NOUT  ,34C)  RHC 

340  FORMAT( •+•  ,83X, 'CORRELATION  COEF F IC I ENT • , 13XF9 . 6 ) 

CALL  OUTPUT ( CELL ( 1  ,  4),  SCALE,  4,  0.) 

WRITE(N0UT,35C)  R 

350  FORMAT ('  +  ' , 83  X  , ' SPEARMAN' •$  RANK  CORRELATION  CCEF. 

C , F9 . 6 ) 

CALL  REPEAT 

IF  (  N  .LE.  5 CO  )  WPITE(NOUT,360)  T 
360  FORMAT( •+',  83X, 'KENDALL"  S  TAU  C C E FF I C I  ENT • , 1 IX F9 .6  ) 
CALL  0UTPUT(CELL( 1  ,  5),  SCALE,  5,  C.) 

CALL  REPEAT 

CALL  OUTPUT ( C  ELL ( 1 ,  6),  SCALE,  6,  Y  M A  X  -  5 .  #  DELTA) 
WR I T  E ( NOUT ,370 ) 

370  FORMAT! •  +  • .95X, • TESTS  FOR  FQU I D IS T P  I8UT ION  '  ) 

CALL  REPEAT 

CALL  OUTPUT (C  ELL ( 1 ,  7),  SCALE,  7,  0.) 

WRITE!  NOUT  ,3  80 

380  FORMAT! '  +  ', 11 IX,  'TEST'  ,  7X, *  NORMALIZED* ) 

CALL  REPEAT 
WR  I  TE  (  NOUT  ,390 

390  FGRMAT( '  +  ' , 109X, *  STATISTIC'  , 5 X , • ST  AT  I  ST  I C ’ ) 

CALL  OUTPUT ( C ELL ( 1  ,  8),  SCALE,  8,  0.) 

WiR  I  TE  (NOU  T,  40  C  )  D 

400  FORMAT! •  +  •  ,83X,» KOLMCGC FOV-SM I RNOV  TEST  »,F10.6) 

CALL  REPEAT 

WR  I  TE ( NOU  T , 41 C )  RUNS,  RNORM 

410  FORMAT! '+ • , 83X  ,  • WALD-WCLFOWI TZ  RUNS  TEST '  ,  1 10 , 5XF9 . 5  ) 
CALL  OUTP  UT ( C  ELL ( 1 ,  9),  SCALE,  9,  0.) 

WRITE! NOUT ,42C)  U,  UNORM 
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420  FQRMAT<  •  +  •  ,  83  X  ,  ’  MIA  NN-WH I TNEY  TEST •  ,  7X  110, 5XF9  .5  ) 


430 


440 


450 


460 


470 


480 


490 


500 


510 


520 


520 


CALL  REPEAT 
WRITE!N0UT,43O 


FORMAT! •+  * ,83X 


CELTA) 
,7X110,  5XF9 .5 ) 

0 


WN,  WNCRM 

WILCOXSON  TEST' ,  10X1 10 ,5X>=9 . 5  ) 
CALL  OUTPUT  ( C  ELL  ( 1  ,  10 )  *  SCALE,  10,  YM  AX  -  9.  * 
WRITE! NOUT  ,  44 C  )  SN,  SNCPM 
FORMAT! •  +  ’ ,  83  X, ' ST  E  GEL— TUKEY  TEST ' 

CALL  REPEAT 

CALL  OUTPUT(CELL( 1  ,11) ,  SCALE, 11,  0.) 

CALL  REPEAT 
WRITE! NOUT  ,  45 C  ) 

FORMAT! '  +  •  ,92  X , ' UN! VAR! ATE  STATISTICS’  ) 

CALL  OUTPUT(CELL( 1,12),  SCALE, 12,  C.) 

CALL  REPEAT 
WRITE!  NOU  T  ,  46  C ) 

FCPMAT! •+• ,95X, ’X  SAMPLE' ,9X, ' Y  SAMPLE') 

CALL  OUTPUT(CELL( 1 ,13) ,  SCALE, 13,  0.) 

WRITE! NOU  7,470  XMOM(l),  YMOM(l) 

FORMAT! •+•  , 83 X, 'MEAN', 5 X1PE14 .6 , 3 X E 14 .6 ) 

CALL  REPEAT 

WRITE! NOU T , 48 C )  XMEC,  Y MED 

F0PNAT(*+',83X,'MEDIAN',3X1PE14.6,3XE14.6) 

CALL  OUTPUT (CELL! 1,14),  SCALE,  14,  Y M A X  -13. 
CALL  REPEAT 

WRITE! NOUT, 490  XM0M<2),  YM0M(2) 

FORMAT! •+•  ,83X,'VARIANCE',1X1PE14.6,3XE14.6) 
CALL  0UTPUT(CELL(1  ,15),  SCALE, 15,  C.) 

WRITE! NOU T , 50C )  XOEV,  YOEV 

FORMAT! •+ • ,83X, 'STD  DE V • , 2X IP  El  4. 6 , 2 XE14. 6 ) 
CALL  REPEAT 
WRITE!  NOUT  ,510 


*  DELTA) 


tixi  icunuui  ,3iu  XRANGE,  YRANGE 
FCRMAT('+',83X,'RANGE',4X1PE14.6,3XE14.6) 
CALL  OUTPUT(CELL( 1 ,16) ,  SCALE, 16,  0.) 

CALL  REPEAT 

WRITE! NOUT, 520  XMCM(4),  YM0M(4) 

FORMAT! *+*  ,83X,' SKEWNESS*  ,1X1 PE14.6,3XE14. 6) 
CALL  OUTPUT (CELL(1»17),  SCALE, 17,  0.) 

. . J  _  vMHM  I  A  1 

6 ,3XE14.6) 
Y  M  AX  -17. 


CALL  OUTPUT (CELL! 1  , 17) ,  SCALE, 17, 
WRITE! NOUT , 53  C  )  XMCM(6),  YM0M(6) 
FORMAT (•  +  •  ,83X,' KUR70S! S'  ,1X1  PE  14. 


18 


*  CELTA) 


CALL  REPEAT 

CALL  OUTPUT(CELL( 1  ,18) ,  SCALE 
WR  ITE!  NOUT, 540  X(N),  YMAX 
540  FORMAT! •+• , 83 X , • MAX  I MUM ’ , 2X1P El 4 . 6 , 2 XE 14. 6 ) 

CALL  REPEAT 

WRITE (NOUT  ,55C )  X(l),  YMIN 
550  FORMAT! •+• ,  83 X, 'MINIMUM ',2X1PE14.6,2XE14.6) 

Z  =  YMAX  -  17.  *  DELTA 
DC  560  J  =  19,32 
Z  =  Z  -  DELTA 

CALL  OUTPUT(CELL( 1, J),  SCALE,  J,  Z) 

CALL  REPEAT 
560  CONTINUE 

WR  ITE  (NOUT ,310 
XLAEEL(l)  =  X(l)  +  OELTAX 
FQUPO  =  4.  *  CELTAX 
CO  570  1=2,8 

XLABEL! I )  =  X  LABEL ( I— 1 )  +  FCURD 
570  CONTINUE 

WRITE (NOUT, 580 (XLABEL! I ) ♦ I =1 , 7, 2 )  ,  ( XLA8FL ( J  )  ,  J  =  2 , 8 , 2  ) 
580  FORMAT! 16X,8(  » I • ,7X)/  12X  ,4! 1PE10 .3  ,  •  I  •)/  20X,4(5 
CIO. 2 ,6X) ) 

KEY! 1 )  =  0 
CO  590  1=2,16 
KEY  (I  )  =  SCALE  *  T  -  1 
590  CONTINUE 

WRITE (NOUT, 600  <KEY(I )  ,1  =  1,16) 

600  FORMAT ( // 50 X, 'KEY' , //'  SYMBOL  PRINTED', 12X 
1-  =  =  =  =  T  > 

2  F  H  F  F'/  • +’  , 38X  ,  •  . 

3  +  X  <  X  V  E  E 

4/  '  + '  ,  56  X ,  '  . 

5/  T  X'/  •+'  ,11  OX,  'T'// 


•  _ 
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6  •  NO.  OBSERVATIONS' , 15, 1516 ) 

RETURN 

END 

FUNCTION  WEIGHTU,  N ) 

INTEGER  WEIGHT 

THIS  FUNCTION  SUBPROGRAM  COMPUTES  THE  THE  WEIGHT  FUNCT 
ION  FOR  THE  I-TH  RANK  IN  A  S I EGEL-TUKEY  TEST  WITH  A  SA 
-MPLE  OF  N. 

K  =  NOD ( 1,2) 

WEIGHT  =  I  +  I  -  K 

I F ( I  .GT.  N)  WEIGHT  =  4  *  N  -  WEIGHT  +  2  *  (1  -  K) 

RETURN 

END 

SUBROUTINE  OUTPUT  (COUNTS,  SCALE,  L,  Y) 

THIS  SUBROUTINE  OUTPUTS  A  SINGLE  LINE  OF  THE  BIVARIATE 
HISTOGRAM  ON  THE  LINE  PRINTER 
ON  THE  LINE  PRINTER 

INTEGER*2  COUNTS,  SCALE,  ZERO,  ONE 
LOG  ICAL*1  LINE,  CHAR 

DIMENSION  L I NE (64 , A ) ,  CHAR(4,  16),  COUNTS ( 32  )  ,  N(16), 
CCH A ( 16) 

EQUIVALENCE  (CHAR(1,1),  C HA (  1 )  ) 

DATA  ONE/1/,  ZERO/O/,  NOUT/6/,  N/3*  1 , 3* 2 , 2  +  3, 2 , 6 

C*3, 4/ 

CATA  CHA  /'  ','  =  ','=.  ','=1 

1  ',  '=+.  ' , ' T X  ',•><-  ','TX- 

2','AVT  »,'HE=  ','HE/  ',  '  HET  ','HEXT'/ 

LINES  =  l 
K  =  1 

CO  20  J=1  ,32 

INDEX  =  CGUNTS(J)  /  SCALE  +  ONE 

IF ( INDEX  .EG.  1  .AND.  CCUNTS(J)  .GT.  ZERO)  INCEX  =  2 
IF ( N ( INDEX )  .GT.  LINES)  LINES  =  N(INCSX) 

DO  1C  1=1,4 

LINE(K,I)  =  CHAR ( I,  INOEX) 

LI NE ( K+l » I )  =  CHAR (  I  , INDEX) 

10  CONTINUE 
K  =  K  +  2 
20  CONTINUE 

IF  (M0D(L-2,4)  .EQ.  0)  GO  TO  40 
ENTRY  REPEAT 

WRI TE ( NOUT , 30  )  (<LINE(J,I),  J=l,64),  1=1, LINES) 

30  FCRNAT(13X,  ’  |  ' , 64 A1  ,  '  |  ’/(•  +  ', 13X , 64 A 1 )  ) 

GO  TO  60 

40  WRITE (NOU  T, 50  )  Y,  (<LINE(J,I),  J  =  l,64),  1  =  1, LINES) 

50  FOPMATC 1X1PE10.3, '  *',64A1,'*'/  ( ' 4 ' , 13X , 64 A  1 )  ) 

60  RETURN 
END 


SUBROUTINE  MOMENt  (WORK,  XMOM) 

PUR  FOSE  : 

FTNO'thE  MOMENTS  OF  THE  OATA  ARRAY  X  WHEN  NOT  ALL 
OF  THE  ARRAY  IS  AVAILABLE  AT  ANY  GIVEN  TIME. 

USAGE: 

CATA  POINTS  ARE  PASSED  TO  THE  SUBROUTINE  " M"  AT  A  T 
-IMF  WITH  THE  RUNNING  TOTALS  KEPT  IN  THE  DOUBLE  PRECIS 
-ION  WORK  AREA  VECTOR  'WORK*.  PRIOR  TO  THE  FIRST  CALL 
, WORK ( 1 )  SHOULD  BE  SET  TO  0.0.  DATA  IS  THEN  PASSED  BY 
THE  CALL 

CALL  ACCUM ( X ,  M,  WORK) 

AN  INITIAL  MEAN  ESTIMATE  IS  KEPT  IN  WQRK<4)  AFTER  THE 
FIRST  7  DATA  POINTS  HAVE  BEEN  PASSED.  HIGHER  CENTRAL 
MOMENTS  ARE  FOUND  USING  THIS  ESTIMATE.  ONCE  AT  LEAST  7 
X  VALUES  HAVE  BEEN  ACCUMULATED,  RESULTS  MAY  BE  CBTAINE 
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-0  FROM 

CALL  MCMENdWCRK,  XMOM) 

PESLLTS  ARE  RETURNED  AS  FOLLOWS: 

XMOM(l)  THE  SAMPLE  MEAN 

XMOM(2)  THE  SAMFLE  VARIANCE  (UNEIASEO  ESTIMATOR) 
XMOM (3 )  THE  SAMPLE  THIRD  CENTRAL  MOMENT  ( IN  ei AS  ED 
ESTIMATOR) 

XM0M<4)  THE  SAMPLE  COEFFICIENT  CF  SKEWNESS 
X  MOM( 5 )  THE  SAMPLE  FOURTH  CENTRAL  MOMENT  (UNBIASED 
ESTIMATOR ) 

XMCM ( 6  )  THE  SAMPLE  COEFFICIENT  CF  KURTQSIS 


IMPLICIT  REAL*8  (D) 

Rc A  1*8  WCRK 

DIMENSION  WORK ( 8 ) ,  XMCM(6),  X(M) 

C 

IF(WORK(l)  .LT.  7.0D0)  RETURN 

C - FINC  CORRECTICN  FACTORS 

CMEAN  =  <  WORK ( 2 )  +  WCRK(B))  /  WORK(l) 

CC  =  DMEAN  -  WOR K ( 4  ) 

CC2  =  OC  *  CC  *  WORK(L) 

C - DETERMINE  CORRECT  CENTRAL  MOMENTS 

DM2  =  WORK ( 5 )  -  DC2 

CM3  =  WOR  K ( 6 )  +  WORK ( 7  )  -  DC  *  (3. CO  *  WORK(5)  -  DC2  - 
CCC2  ) 

DM4  =  WOR  K  (  8 )  -  DC  *  (4. DO  *  (WORKU)  +  WORK(7))  -  DC 
1  *( 6 . CO  *  WCRK(5)  -  3. DO  *  C C2 ) ) 

C - CORRECT  ESTIMATORS  IN  WORK  ARRAY 

WOR  K ( 4 )  =  CMEAN 
WOR  K ( 5 )  =  DM2 
WOR  K ( 6 )  =  0. 

WCR  K ( 7 )  =  0. 

IF ( CM3  .GT.  0.)  WORK (6 )  =  DM3 
IF  <  CM3  .LT.  0.)  WORK ( 7 )  =  DM3 
WOR  K ( 8 )  =  0M4 

C - RETLRN  STATISTICS 

XMCM(I)  =  OMEAN 
CF  AC  =  WORK ( I  )  -  1 .DO 
XMCM ( 2 )  =  0M2  /  DFAC 
CF  A  C  =  DFAC  *  ( WORK ( 1  )  -  2 . DO ) 

XMCM ( 3 )  =  WORK ( 1 )  *  DM3  /  DFAC 
XMOM ( 4 )  =  XMC  M ( 3 )  /  <XM0M(2)  * 

CF AC  =  DFAC  *  (WOR  K  ( 1 )  -  3. DO) 

XMO M ( 5 )  =  (  ( ( WORK ( 1  )  -  2. DO)  * 

I  -  (6. DO  *  WORK ( I  ) 

2W0RM1))  /DFAC 

XMOM (6)  =  XMOM ( 5)  /  (XMCM(2)  * 

RETURN 
C 
C 

ENTRY  ACCUM ( X  t  M,  WCRK) 

C 

I F ( WORK < I  )  .GT.  6. DO)  GO  TO  210 

C - ACCUMULATE  FIRST  7  CATA  POINTS 

N  =  WORK ( 1)  +  2. DO 
WORK (  1 )  =  WORK(l)  +  M 

NP  =  |s  +  y  _  i 

IF(NP  .GT.  8)  NP  =  8 
J  =  1 

CO  201  I=N,  NP 
WORK (  I)  =  X ( J ) 

J  =  J  +  1 

201  CONTINUE 
IF  ( N  F  .LT.  8)  RETURN 

C - CBTAIN  INITIAL  MEAN  ESTIMATE 

CSUMP  =  0. 

CSUMM  =  0. 

CO  202  I  =  2,8 

I c ( WORK ( I  )  .LT.  0.)  DSUMM  =  DSUMM  +  WORK(!) 

IF ( WORK ( I )  .GT.  0.)  DSUMP  =  DSUMP  +  WORK ( I ) 

202  CONTINUE 


SORT ( XMOM( 2 ) > ) 

WORK ( 1 )  +  3  .n0  )  *  DM4 
-  S.CC)  *  DM2  *  CM2/ 

XMOM ( 2  )  )  -  3.0 
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IF  ( J  .GT.  M)  GO  TO  204 
CO  203  I=J,M 

I F ( X  ( I  )  .LT.  C.)  OSUMM  =  DSUMM  +  X(I1 
IF ( X (  I  )  .GT.  C.)  OSUMP  =  OSUMP  +  X ( I J 
2C3  CONTINUE 

204  CMEAN  =  (CSUMM  +  OSUMP )  /  W0RK(1) 

C-- « SET  up  WORK  AREA 

0  2  =  0. 

C3M  =  0. 

D3P  =  0. 

04  =  0. 

CO  205  I =  2  t  8 
Cl FF  =  WORK ( I  )  -  DMEAN 
CF2  =  DIFF  *  C IFF 
02  =  02  +  DF2 

I- ( Cl FF  .UT.  C.)  D3N  =  C3M  +  DIFF  *  CF2 

IF  ( C  I  FF  .GT.  C.)  D3P  =  C3P  +  DIFF  *  DF2 

04  =  04  +  Dc2  *  0 P 2 

205  CONTINUE 

IF ( J  .GT.  M )  GO  TO  207 
00  206  I  =  J  ,  M 
CIFF  =  X( I)  -  CMEAN 
Cp 2  =  DIFF  *  CIFF 
02  =  02  +  DF2 

I F  <  C I FF  .LT.  0.)  03  m  =  C3M  +  DIFF  *  CF2 

IF  CC I  FF  .GT.  C.)  03  P  *  C3P  +  DIFF  *  CF2 

04  =  04  +  DF2  *  0C2 


206  CONTINUE 


207  WORK ( 2 ) 

= 

DSUMP 

WORK (3) 

= 

CSUMM 

WOP  K(4) 

55 

DMEAN 

WORK (5) 

55 

D  2 

WCR  K ( 6 ) 

= 

03P 

WORK  ( 7) 

55 

D3M 

WORK (8) 

55 

D4 

RETURN 

C - accumulate  DATA 

210  WCRK(l)  =  WOPK(l)  +  M 
CO  211  1  =  1,  M 

I F ( X (  I  )  .GT.  C.)  WORK ( 2  )  =  W0RK(2)  ♦  X ( I ) 

I F  (  X (  I  )  .UT.  0.)  WORK( 3  )  =  W0RK(3)  +  X  < I ) 

CIFF  =  X(  I)  -  WOR  K ( 4  ) 

CF2  =  DIFF  *  CIFF 

WORK ( 5 )  =  WOP  K ( 5 )  +  DF? 

IF  ( 0 1  FF  .GT.  C.)  WORK  ( 6 )  =  W0RK<6)  4  OIFF  *  DF2 

IF(OIFF  .LT.  C.)  WORK ( 7 )  =  W0RK(7)  +  OIFF  *  DF2 

WOR K ( 8 )  =  WORK (8 )  +  CF2  *  0F2 

211  CONTINUE 
RETURN 
END 

C 

C 

SUBROUTINE  SCRTON  (ON,  II,  JJ,  WITH 

THIS  SUBROUTINE  SORTS  TFE  VECTOR  ON  INTO 
ASCENDING  ORDER  FROM  CN  (II)  TO  ON(JJ),  CARRYING 
ALONG  THE  CORRESPOND  ING  ELEMENTS  OF  THE  VECTOR 
WITH 

INTEGER  HI,  HISTK,  FIKEEP,  ST  KPT ,  WITH 
DIMENSION  ON(JJ),  WITH(JJ),  HISTK(16),  L0STK(16) 

STKPT  =  1 
LOKEEP  =  II 
FIKEEP  =  JJ 

10  IF  (LOKEEP  .GE.  HlKEEP)  GO  TO  100 
20  LO  =  LOKEEP  -  1 
FI  =  HlKEEP  +  1 

MIDCLE  =  (LOKEEP  +  HlKEEP)  /  2 
IL  =  LOKEEP 
IH  =  HlKEEP 

IF  (  ON ( MIDDLE )  .GT.  ON (HIK  SEP )  )  IF  =  MIDDLE 
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IF 

( 

ON ( MIOOL  E  ) 

.LT.  ON ( LOKEEP ) 

) 

IL  =  MIDOLE 

IF 

( 

0N( HIKEEP) 

.LT.  ON(IL) 

) 

IL  =  HIKEEP 

IF 

( 

ON ( LOK  EE P ) 

.GT.  ON(IH) 

) 

IF  =  LOKEEP 

ISUE 

=  LOKEEP  + 

MIDCLE  +  HIKEEP 

— 

IF  -  IL 

TEST  =  ON(ISUe) 
GO  T0  50 


30  IF  (LO  .EQ.  HI)  GO  TO  50 
TEMP  =  CN(HI) 

ON ( F I  )  =  CN(LC) 

CN(LO)  =  TEMP 
ITEMP  =  WI  TH(  HI  ) 

MTF(HI)  =  WITH(LO) 

WITF(LO)  =  ITEMP 
50  HI  =  HI  -  1 

IF  (  ON ( H I )  .GT.  TEST)  GO  TO  50 
70  LC  =  LO  +  1 

IF  (  CNILO)  .LT.  TEST)  GO  TO  70 

IF  (  LO  .LE.  HI  )  GO  TO  30 

I F ( F I  -  LCKEEP  .LE.  HIKEEP  -  LO)  GO  TO  90 

LOSIK(STKPT)  =  LOKEEP 

FISTK(STKPT)  =  HI 

STKPT  =  STKPT  +  1 

LOKEEP  =  LO 

GO  TC  110 

90  LOSTK(STKPT)  =  LO 

HI STK ( STKPT )  =  HIKEEP 
STKPT  =  STKPT  +  1 
FIKEEP  =  FI 
GO  TO  110 

100  STKPT  =  STKPT  -  1 

IF  < STKPT  ,£Q.  0)  RETURN 
LCKEEP  =  LOSTK(STKPT) 

HIKEEP  =  HI  ST  K (STKPT ) 

110  IF  (  HIKEEP  -  LOKEEP  .GE.  11)  GO  TC  20 
IF  (  LOKEEP  .EQ.  II)  GO  TO  10 
LOKEEP  =  LOKEEP  -  1 
120  LCKEEP  =  LCKEEP  +  1 

IF  (  LOKEEP  .EQ.  HIKEEP  )  GO  TO  ICO 

IF  <  ON(LCKEEP)  .LE.  CN(L0KESP+1)  )  GO  TO  120 

TEMP  =  ON ( LOK  EEP  +  1) 

LQ  =  LCKEEP 
130  CN (  10+ 1 )  =  CN(LO) 

LC  =  LO  -  1 

IF  (  TEMP  .LT.  CN(LC))  GO  TO  130 
CN ( LO+1 )  =  TEMP 
I VI  =  LCKEEP  +  1 
N  =  LCKEEP  -  LO 
ITEMP  =  WITHIIW) 

CO  140  IMOVE  =  1 »  N 
hlTF(lW)  =  WITH(IW-l) 

IW  ■  IW  -  1 
14C  CONTINUE 

WITF(IW)  =  ITEMP 

GO  TO  120 

ENO 
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onnn  o  r^ooooor*Oooor-YooooooOoooo<~> 


FCP  GENERATING  EIVARIATE  RANDOM  V  ECTOR  S  (  Y,  Z  )  ♦  C  A  LL 
FOLLOWING  SUBROLTINE . 

1.  UN  I  FOR  M  (  0 , 1 )  MARGINAL  DISTRIBUTION. 

CALL  LEWUFM(CR,NV.Y,Z,IS);  POSITIVE  CORRELATION  BY 
LAWERENCE  AND  LEWIS  METHOD. 

CALL  GAVUFM(CR  ,NV,Y,  Z,  IS  )  ;  NEGATIVE  CORRELATION  BY 
TRANSFORMATION  OF  GAVER’S  EXPONENTIAL. 

2.  EXPONENT  IAL(MEAN=I.Q)  MARGINAL  Cl  STRIBUTICN  • 

CALL  MAREXF(CR,NV,Y,Z,IS) ;  POSITIVE  CORRELATION  BY 
MARSHAL  ANC  OLKINS  METHOD. 

CALL  GAVEXF(CR,NV,Y,Z, IS);  NEGATIVE  CORRELATION  BY 
GAVER'S  METHOD. 

WHERE 

CR;  GIVEN  CORRELATION. 

NV;  NUMBER  OF  RANDOM  VECTORS  BE  GENERATING. 

Y;  THE  RIRST  ELEMENT  OF  VECTOR  WITH  CIMENSICN  NV. 

2;  THE  SECOND  ELEME  NT  OF  VECTOR  WITH  DIMENSION  NV. 

IS;  INITIAL  RANDOM  SEED. 


SUBROUTINE  L EWUFM ( CR ,NV ,Y , Z , IU ) 

DIMENSION  Y (NV ) ,Z ( NV ) , U ( 3) 

IF  (CR.LT  .0.0)  GO  TC  20 

A  =  S  CRT (36.0+(  12. 0*  CR  )  +  (  16.0*(CR**2  .C  )  ))-6  .0 
6 ET A  =  A/ ( 2 . 0*C P  ) 

ALP A=2.0-(I.0/8ETA ) 

F=(  l.C-BETA) /(1.0-BSTA+(ALPA=8ETA ) ) 

R  =  1  .O-ALP A 

CO  ICO  1=1, NV 

CALL  RANDOM ( IU,U,2) 

Y  (  I  )  =  U(1  ) 

IF ( L ( 2) .LE.R)  GO  TC  50 
W= (U(2)-R )/ ( I.O-R ) 

IF(W.LE.P  )  GO  TO  10 

2(I)=(Y(I)**BETA)*<<W-P)/(i.0-P))**(P*eETA) 
C-C  TC  100 

10  Z(I  )  =  (Y(I)**eETA)*(W/P ) 

GO  fC  100 
5 C  V  =  U  (2  )/R 

IF(  V.LE.P )  GO  TO  60 

Z(  I)  =  ((W-P)/( 1.0-P) ) ** ( R*BE  TA ) 

GO  TC  100 
6  C  Z(I)=W/P 
100  CONTINUE 
PETLRN 

2C  WRITE(6,30)CR 

30  FCRMAT( 5X, 'THIS  METFCD  NOT  ALLOW  CR=',F7.3) 
PETLRN 
ENC 


SUBROUTINE  GA  VUFM ( C  P , NV  ,Y , Z , I  S) 

C 

DIMENSION  Y(NV),Z(NV),WK(1),G(100) 

COUELE  PRECISION  DXS 

CXS=IS*0.8 

Il=IS*0. 5+12345 

IF  (CP  .GT  .0.0  )  GO  TO  20 

P=-2.0*CR 

C=1 .C-P 

CC  100  1=1, NV 

CALL  GGEOT(CXS ,1 ,0  ,WK1  ,  IG) 

CALL  RANDOM ( II, U, 1 ) 

RN=1 .0/1 G 

Y  ( I  )  =  ( 0  =*  (  N  )  )/<1.0-(P*(U**RN)  )  ) 
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CALL  RANDOM ( It ,G, I G  ) 

V*  =  1 . 0 

CG  1C  K=1 , I G 
W=W*G<  K  ) 

1C  CGMINLE 

Z< I )=W**( l.O-F) 

ICO  CONTINUE 
R  ETLRN 

20  WRIT5<6,30)CP 

30  FCRNAKSXi’THlS  METhCO  NOT  ALLOW  CR  =  '  » F7 . ?  ) 
PETLPN 
END 
C 
C 
C 
C 

SUBROUTINE  M A R EXP < CR , NV , Y, Z , I  X ) 

C 

DIMENSION  Y(NV),Z(NV),E(3) 

Ip  (CR.LT.O.O)  GO  TC  20 
R12X2.C*CR)/(CR+1.C) 

P 1=  1 .  0-R 1 2 
P2  =  R  1 

DO  5  C  1=1, NV 
CALL  EX  PON  ( I  X  ,  E,  3  ) 

E (1 )  =  E< 1) /R1 
E<  2  )  =  E(  2)  /R2 
E  (3  )*E(3)/R12 
X I ) =AM I N 1 ( E  < 1 )  ,6(3)) 

Z(  I ) =  AM IN 1( E( 2) ,E( 3  )  ) 

50  CONTINUE 
RETURN 

20  WRITE(6,3C)CR 

30  FCRVAT( 5X  ,  • ThIS  METhCD  NOT  ALLOW  CR=’,F7.3) 
PEURN 
END 


C 

C 

SUE  ROUTINE  C-AVEXP  (  CR  ,  NV  ,  Y,  Z  ,  I  S  ) 

C 

DIMENSION  Y(NV) ,Z( NV) , WK1( 1) , WK2( ICC) 

CCLELE  PRECISION  DS EED1 ,0SEED2, DS EEC3 
CSEED1=IS*0. 2  +  123456 
0SEEC2=IS*0. 5+1234 
CSEED3=IS*0  .9  +  123 

IFICR.LT. -0. 5  .OR.  CR.GT.O.O)  GO  TC  20 
C  =  -2  .0*CR 
P*1.C-Q 
DO  ICC  1=1, NV 

CALL  GGEOT< CSEED1 ,1,P,WK1, IG) 

A=  I G 

e=p 

CALL  GGAMS(CSEED2,A,B,1,WK2,GA) 

Z(I ) =GA 

CALL  GGUBS( DSEED3, 1  ,U) 

PN  =  1 .0/  IG 

XI  )=ALOG  ((  1.  C/(  »*U**RN)  )-<  C/F)  ) 

100  CONTINUE 
RETLRN 

20  VRITE(6,3C)CR 

30  FGPWAT(5X,»THIS  METhCD  NOT  ALLOW  CR=',F7.3) 
RETURN 
END 
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